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0.1  Abstract 


This  technical  note  investigates  the  estimation  of  environmental  parameters  in  the 
Bistatic  Scattering  Strength  Model  (BISSM)  given  backscatter  strength  and  bathy¬ 
metric  data.  A  monostatic  version  of  the  model  is  derived,  since  this  will  be  the  form 
of  data  provided  by  acoustic  imaging  sensors.  Feedforward  neural  networks,  using  the 
backpropagation  learning  algorithm,  are  used  to  perform  the  estimation  of  parameters 
for  the  nonlinear  BISSM  equation.  The  parameters  that  can  be  estimated  are  iden¬ 
tified,  and  neural  networks  have  been  developed  to  estimate  these  parameters.  Using 
noise-free  artificial  data  generated  with  the  BISSM  equation,  the  networks  provided 
excellent  estimates  of  the  desired  parameters. 

The  primary  impetus  for  this  work  is  a  need  for  the  Naval  Oceanographic  Office 
(NAVOCEANO)  to  provide  relevant  survey  support  for  Low-Frequency  Active  Acous¬ 
tics  (LFAA)  programs  and  future  Low-Frequency  Active  (LFA)  operational  systems. 
It  has  been  recognized  by  the  Commander,  Naval  Oceanographic  Command  (CNOC) 
that  such  support  will  require  knowledge  of  certain  bottom  and  subbottom  properties 
and  high-resolution  geomorphology. 

The  BISSM  model  has  been  proposed  as  a  model  for  aspects  of  active  bottom  rever¬ 
beration.  It  was  postulated  that  the  parameters  that  activate  the  BISSM  algorithm, 
might  be  measurable  with  NAVOCEANO’s  swath  bathymetry  system  -  SASS  phase 
IV.  This  latest  version  of  the  SASS  system  developed  by  the  Naval  Air  Development 
Center  (NADC)  generates  91  one-degree  beams  that  can  reach  grazing  angles  down  to 
45  degrees  and  records  backscattering  strength  as  a  function  of  grazing  angle.  Future 
SASS  systems  are  expected  to  reach  grazing  angles  down  to  30  degrees.  Although 
BISSM  is  designed  as  a  LFAA  prediction  tool  and  should  be  most  valuable  as  such 
below  IkHz,  it  should  be  scaleable  to  the  higher  frequency  of  12kHz  used  by  SASS 
and  useful  for  inverse  applications.  The  system  parameters  of  S.4SS  will  determine 
the  characteristics  of  the  measurement  process. 
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Neural  Network  Parameter  Estimation  for  the  Modified 
Bistatic  Scattering  Strength  Model  (BISSM) 


Chapter  1 
Introduction 


This  paper  investigates  estimation  of  the  environmental  parameters  in  the  Bistatic 
Scattering  Strength  Model  (BISSM).  The  BISSM  model  has  recently  been  proposed  by 
Caruthers  et  al.  [1]  as  an  advancement  of  the  model  previously  developed  at  the  Naval 
Oceanographic  and  Atmospheric  Research  Laboratory  (NOARL).  The  end  goal  is  to 
use  the  backscatter  and  bathymetric  data  obtained  from  advanced  acoustic  imaging 
sensors  to  determine  the  parameters  in  the  BISSM  model.  Thus,  given  backscatter 
strength  and  the  incident  and  azimuth  angles  determined  from  bathymetry,  it  is 
desired  to  compute  the  ratio  of  sound  speeds  and  densities  at  the  bottom  interface, 
the  Mackenzie  coefficient,  the  root  mean  square  (RMS)  microscale  heights  roughness, 
and  the  fine-scale  RMS  slopes  in  the  along  track  and  across  track  directions. 

The  approach  taken  in  this  technical  note  is  to  use  artificial  data,  generated  by  the 
BISSM  equation,  to  estimate  the  known  parameters  that  were  used  to  generate  the 
data.  Since  BISSM  is  a  nonlinear  equation  it  cannot  be  inverted  directly,  requiring  an 
error  minimization  approach  for  the  estimation  of  its  parameters.  Using  artificial  data 
provides  an  opportunity  to  determine  the  best  possible  performance  of  the  estimation 
techniques.  Neural  networks  are  used  in  this  work  to  estimate  the  desired  parameters, 
given  an  input  data  vector  of  backscatter  strength  as  a  function  of  incident  angle. 
Essentially,  this  work  shows  proof  of  concept,  using  noise  free  data. 

In  the  next  chapter,  the  BISSM  equation  is  presented,  and  the  monostatic  ver¬ 
sion  of  this  equation  is  derived.  A  monostatic  version  is  used  in  this  work  since 
this  is  the  form  of  data  that  will  be  collected  by  acoustic  imaging  systems.  Chap¬ 
ter  3  investigates  the  sensitivity  of  the  backscatter  strength  produced  by  the  BISSM 
equation  to  its  parameters.  By  examining  the  correlation  between  the  effects  of  the 
various  parameters,  a  set  of  parameters  that  have  potential  for  estimation  are  deter¬ 
mined.  In  Chapter  4  various  multidimensional,  nonlinear  minimization  approaches 
are  discussed,  and  the  advantage  of  a  neural  network  approach  over  more  traditional 
approaches  is  explained.  A  quick  review  of  feedforward  neural  networks  and  the  back- 
propagation  learning  algorithm  is  also  given  in  this  chapter.  Chapter  5  presents  the 
results  of  the  parameter  estimation  attempts,  showing  that  the  trained  neural  net¬ 
works  reliably  provide  good  estimates.  Finally,  Chapter  6  summarizes  the  significant 
results,  and  details  further  research  that  is  required. 
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Chapter  2 

Monostatic  Reduction  of  BISSM 


As  stated  in  the  previous  chapter,  the  goal  of  this  task  is  to  estimate  the  environ¬ 
mental  parameters  in  the  BISSM  equation  given  collocated  bathymetric  and  backscat- 
ter  data.  Specifically,  it  is  desired  to  estimate  n,  the  ratio  of  sound  speeds,  m,  the 
ratio  of  densities,  /i,  the  Mackenzie  coefficient,  cr,  the  RMS  microscale  heights  rough¬ 
ness,  and  <5x  and  6y,  the  fine-scale  RMS  slopes  in  the  x  and  y  directions.  Potential 
sources  of  data  for  this  inversion  are  multibeam  hull  mounted  sonar  and  towed  sides- 
can  sonar  systems.  In  both  cases  the  source  and  receiver  will  be  at  the  same  position 
so  that  6  ~  6i  =  6s  and  4>  =  4>i  =  <i>s  ~  i  where  is  the  source  incident  angle.  Og 
is  the  receiver  scattered  angle,  and  (t>,  and  4>s  are  the  incident  and  scattered  azimuth 
angles.  Thus,  the  first  step  in  parameter  estimation  is  to  redu''e  the  bistatic  model 
to  a  monostatic  model. 

The  BISSM  bistatic  equation  [1]  is  given  by: 

mbs  =  mi  +  m2  .  (2.1) 

In  (2.1)  mi  is  the  incoherent  term  and  is  given  by; 

mi  = /i sin  0,  sin  0s  .  (2.2) 

Both  angles  are  measured  upward  from  the  local  bottom  facet  plane.  In  (2.2)  p  is 
the  Mackenzie  coefficient. 

In  (2.1)  m2  is  the  coherent  term  and  is  given  by; 

m,  =  /!jexp(-,)^exp|-i[^  +  :h)|  .  (2.3) 

Ro,  in  (2.3),  is  the  Rayleigh  reflection  coefficient  between  two  fluids  which  is  given 

by: 

m  sin  0,  —  \/n^  —  cos'^  0, 

Ro  =  - ; — - /  (-•‘1) 

m  sin  Pi  -f  \/n‘  —  cos^  0, 

where  m  is  the  ratio  of  densities  and  n  is  the  ratio  of  sound  speeds.  In  (2.3)  g  is  the 
square  of  the  Rayleigh  roughness  parameter  given  by; 

9  =  (■2-5) 
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where  a  is  the  RMS  microscale  heights  roughness,  q  is  given  by: 


q  =  k{sm  0,  +  sin  6s)  (2.6) 

and  k  is  the  acoustic  wavenumber  given  by  27r/A,  where  A  is  the  wavelength  of  the 
acoustic  source.  F  is  given  by: 

where 

=  xl  +  xl 

and 


Xx  =  k  {cos  6 s  cos  (l>s  —  cos  9 i  cos  <f>i)  (2.8) 

Xy  =  A:  (cos  sin  ^5  —  cos  0,  sin  . 

4>i  and  (^5  are  measured  clockwise  from  North  to  the  projection  of  the  vector  onto  the 
local  horizontal  plane.  Finally,  Sx  and  Sy  are  the  fine-scale  RMS  slopes  in  the  x  and 
y  directions,  and  are  given  by: 


inn  « 

”  .=lj=l 

(2.9) 

inn  « 

.=ij=i 

(2.10) 

where  n  is  the  number  of  pixels  along  one  side  of  a  square  grid,  and  6'J  and  are 

the  incremental  change  in  slope  in  the  x  and  y  directions.  The  equation’s  parameters 
are  summarized  in  Table  2.1. 

For  the  monostatic  Ccise,  the  incoherent  term  (2.2)  simplifies  to: 

mi  =  fi  sin  9  sin  9 

=  /isin^  9  . 

(2.11) 

Also,  Xx  becomes 

Xx  =  A:  [cos  0  (cos  <?i>  —  cos  (<^  —  7r))] 

=  k  [cos  9  (cos  4>  -f  cos  (j))] 

=  2k  (cos  9  cos  <f)) 

and  Xy  becomes 

Xy  =  k  [cos  9  {sin  <j>  —  sin  {(f)  —  tt))] 

=  A:  [cos  (sin  ^ -t- sin  <?!>)] 

=  2k  {cos  9  sin  <f>)  . 
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Table  2.1:  Parameters  in  the  BISSM  model. 


Parameter 

Description 

Units 

n 

ratio  of  sound  speeds 

unitless 

m 

ratio  of  densities 

unitless 

fi  (mu) 

Mackenzie  coefficient 

unitless 

a  (sigma) 

RMS  micro-scale  heights  roughness 

meters 

Sx  (delx) 

fine-scale  RMS  slopes  in  x  direction 

radians 

Sy  (dely) 

fine-scale  RMS  slopes  in  y  direction 

radians 

9. 

source  incident  angle 

radians 

9s 

receiver  scattered  angle 

radians 

9 

monostatic  incident  angle 

radians 

4>i 

source  azimuth  angle 

radians 

4>s 

receiver  azimuth  angle 

radians 

k 

wavenumber 

radians/meter 

q  reduces  to 

and  from  (2.3)  we  have: 


q  =  2k  sin  d 


(  1 


=  exp  < - 


1  f  4k^  cos^  0  cos^  (f)  ^  Ak^  cos^  6  sin^  q'' 


8k'^  sin^  9 


SI 


f  cot^  0  /  cos^  </>  sin^ 
=  exp  < - ^ —  1  —t::: - 1- 


V  2  V 

Finally.  F  can  be  simplified  to 

F  =  + 


H 


1  /  cos^  0  cos^  0  +  cos^  0  sin^  <A' 


2"  + 


Ak"^  sin^  6 


\  /  4P  cos^  0  ^cos^  0  +  sin^  (i^j  \ 
2  ^  ^  Ak'^sin^O  j 


- 

=  ^  csc^  6 


cos^  0 
sin^  9 


and 


= 


1 


4  sin'*  9 


A 


Thus,  (2.1)  becomes 


Rlexp{-g) 


cot^  9 


cos^  (j)  sin^  <f) 


Note  that  if  we  assume  6  =  =  Sy  then  (2.12)  simplifies  to 


exp 


(- 


cot^  9 

2P 


and  (2.13)  simplifies  to 

TUma  =  /xsin''^  + 


2z.  ,  ^exp(-</) 


sin^  9 


exp 


cot^  9 
■  2<52 


(2.13) 


(2.14) 


Given  collocated  bathymetric  and  backscatter  data  for  a  region,  we  desire  to 
determine  the  values  of  n,  m,  /i,  a,  and  8y  that  provide  the  least  mean  square  error 
between  the  actual  data  (for  a  single  scan  line)  and  (2.13).  The  angles  9  and  (f)  are 
determined  from  the  bathymetric  data  and  the  source  ray  trace  angle  (with  respect 
to  the  local  horizontal  plane),  and  k  is  determined  from  the  source  frequency  and  the 
local  speed  of  sound  in  the  water  (assumed  constant). 


Chapter  3 


Independence  of  BISSM 
Parameters 

As  previously  staled,  it  is  desired  to  estimate  the  parameters  n,  m,  /x,  cr,  and 
Sy  in  (2.13)  given  collocated  bathymetric  and  backscatter  data.  Successful  estimation 
of  these  parameters  requires  that  their  effect  on  the  computed  backscatter  intensity 
be  uncorrelated  with  each  other.  In  this  chapter  the  correlation  between  the  desired 
parameters  is  determined  empirically  in  order  to  identify  those  parameters  that  may 
potentially  be  estimated.  The  ranges  of  interest  for  the  desired  parameters  for  abyssal 
plains  are  given  in  Table  3.1.  A  wavenumber  of  5.0265  radians/meter  is  used,  which 
corresponds  to  a  sound  velocity  of  1500  m/sec  and  a  sonar  frequency  of  1.2  kHz. 


T  able  3.1:  Parameter  ranges  for  the  BISSM  mod  el 


Parameter 

Nominal 

Low 

High 

n  (unitless) 

0.99 

0.97 

1.2 

m  (unitless) 

1.4 

1.2 

1.8 

/X  (unitless) 

0.002 

0.0002 

0.02 

cr  (meters) 

0.01 

0.005 

0.02 

Sx  (radians) 

0.08727 

Sy  (radians) 

0.05236 

Figure  3.1  shows  the  backscatter  strength  versus  incident  angle  with  all  six  param¬ 
eters  at  their  low,  nominal,  and  high  values.  In  this  figure  the  angle  has  been  set  to 
zero.  Figure  3.1  illustrates  three  dominant  motions  in  the  backscatter  strength  curve 
as  the  parameter  values  are  varied.  One  motion  is  a  fairly  uniform  rise  in  backscat¬ 
ter  intensity  at  lower  incident  angles  (less  than  70  degrees)  as  the  parameters  are 
increased.  Secondly,  the  ‘pivot  point’,  where  the  backscatter  intensity  curves  sharply 
upward,  moves  to  lower  incident  angles  as  the  parameter  values  are  increased.  Thirdly, 
the  slope  of  the  sharp  rise  at  higher  angles  decreases  with  an  increase  in  parameter 
values. 
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Incident  Angle 


Figure  3.1:  Scattering  Strength  Curves  for  all  parameters  at  their  high,  nominal,  and 
low  values. 


Figure  3.2:  Scattering  Strength  Curves  for  //  at  high,  nominal,  and  low  values,  with 
all  other  parameters  nominal. 


0 


mu  =  0.002 


Figure  3.3:  Scattering  Strength  Curves  for  n  at  liigh,  nominal,  and  low  values,  with 
all  other  parameters  nominal. 

Figure  3.2  shows  the  effect  on  backscatter  strength  due  to  variations  in  with  all 
other  parameters  at  their  nominal  values.  It  is  seen  in  this  figure  that  an  increase  in 
l-i  will  result  in  a  fairly  uniform  increase  of  the  backscatter  strength  for  lower  incident 
angles  (less  than  70  degrees).  From  the  remaining  figures  shown  in  this  chapter,  it 
is  seen  that  is  the  dominant  parameter  for  backscaHer  intensity  changes  at  low 
incident  angles. 

Figures  3.3  through  3.5  show  the  effect  on  backscatter  strength  due  to  variations 
in  the  n  and  in  parameters.  In  Figure  3.3  n  is  varied  and  all  other  parameters  have 
nominal  values,  in  Figure  3.4  m  is  varied,  and  in  Figure  3.5  n  and  m  are  varied 
simultaneously.  Both  parameters  cause  a  change  in  the  slope  of  the  sharp  rise  at  high 
incident  angles,  where  an  increase  in  n  causes  a  decrease  in  slope,  and  an  increase 
in  m  causes  an  increase  in  slope.  It  is  seen  from  these  figures  that  the  effect  of 
changes  in  n  and  m  on  the  backscatter  strength  are  highly  correlated,  so  n  and 
rn  cannot  be  independently  estimated.  However,  when  n  is  large  m  is  also  large, 
so  a  simplifying  assumption  will  be  used  during  parameter  estimation  that  the  two 
parameters  are  linearly  related.  This  essentially  reduces  the  parameter  estimation  to 
a  single  parameter,  m  or  n.  and  m  will  be  estimated  in  this  work.  Figure  3.6  shows 
that  the  effect  on  backscatter  strength  due  to  changes  in  cr  are  highly  correlated  with 
those  due  to  changes  in  n  and  m.  Furthermore,  the  variations  due  to  changes  in  cr  are 
extremely  small.  Consequently,  cr  cannot  be  successfully  estimated  for  a  frequency  of 
1.2  kHz,  especially  since  actual  data  will  be  contaminated  by  noise. 

Figures  3.7  through  3.11  show  the  effect  on  backscatter  strength  due  to  variations 
in  Sy,  and  0.  In  Figures  3.7  and  3.8  6  is  zero,  and  it  is  seen  that  Sj-  shifts  the  ]nvot 
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Figure  3.4.  Scattering  Strength  Curves  for  m  at  high,  nominal,  and  low  values,  with 
all  other  parameters  nominal. 


Figure  3.5;  Scattering  Strength  Curves  for  both  n  and  m  at  high,  nominal,  and  low 
values,  with  all  other  parameters  nominal. 


point  and  affects  the  slope  at  high  incident  angles,  while  only  affects  the  slope.  This 
result  is  expected  from  the  form  of  (2.13)  since  sin  <}>  is  zero  when  Q  is  zero.  In  Figures 
3.9  and  3.10  0  is  45  degrees,  so  the  effect  on  backscatter  strength  due  to  changes  in 
Sj;  and  6y  are  equal.  These  figures  shov,-  that  the  effect  of  changes  in  Sr  and  Sy  are 
highly  correlated,  so  that  these  parameters  cannot  be  independently  estimated.  To 
simplify  the  problem,  it  will  be  assumed  that  8  =  6x  =  Sy,  and  this  single  parameter 
will  be  estimated.  As  shown  in  (2.14),  when  we  set  Sr  =  Sy  the  parameter  (p  vanishes. 
In  Figure  3.11  Sr  and  Sy  are  varied  simultaneously.  The  effect  on  backscatter  strength 
due  to  the  S  parameter  is  seen  to  be  a  shift  in  the  incident  angle  of  the  pivot  point, 
and  a  change  in  the  slope  at  higher  angles.  While  there  is  some  correlation  between 
the  S  and  the  m  parameter  as  seen  by  the  change  in  high  angle  slope,  the  change  in 
pivot  point  due  to  S  may  provide  enough  unique  information  to  distinguish  between 
changes  in  backscatter  strength  due  to  S  versus  nm. 

Summarizing,  it  has  been  observed  that  only  three  parameters  may  potentially 
be  estimated  from  actual  backscatter  data  using  the  BISSAI  (monostatic  version) 
relationship:  /r,  m,  and  S.  The  changes  in  the  /j.  parameter  result  in  a  shift  up  or 
down  of  the  backscatter  strength  at  low  angles  of  incidence  (less  t’  an  70  degrees). 
The  effect  of  the  n  and  m  parameters  on  backscatter  strength  were  seen  to  be  highly 
correlated,  so  only  one  of  these  parameters  can  be  estimated.  A  linear  relationship 
will  be  assumed  between  n  and  m  to  perform  the  estimation,  so  only  ??!  will  be 
estimated.  The  primary  effect  of  the  n  and  m  parameters  is  to  change  the  slo])e  of 
the  curve  at  high  incident  angles.  It  was  observed  that  the  effect  of  a  on  backscatter 
strength  was  too  small  to  be  reliably  estimated  at  the  frequenct’  used  in  this  .stud}’. 
Also,  the  effects  of  cr  and  nm  are  highly  correlated.  The  effects  of  Sr  and  Sy  were  seen 
to  be  highly  correlated,  depending  upon  the  value  of  <p.  Consequently,  these  terms 
are  combined  into  a  single  parameter  5.  The  S  parameter  affects  the  slope  of  the 
backscatter  strength  curve  at  high  incident  angles,  which  is  correlated  with  the  7nn 
parameter,  but  it  also  shifts  the  pivot  point  in  the  curve  transition  from  low  to  high 
incident  angles. 
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Figure  3.6:  Scattering  Strength  Curves  for  a  at  high,  nominal,  and  low  values,  with 
all  other  parameters  nominal. 
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Figure  3.7:  Scattering  Strength  Curves  for  at  high,  nominal,  and  low  values,  with 
all  other  parameters  nominal,  and  <^  =  0  degrees. 
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Figure  3.8:  Scattering  Strength  Curves  for  6y  at  high,  nominal,  and  low  values,  with 
all  other  parameters  nominal,  and  <t>  =  0  degrees. 


Figure  3.9;  Scattering  Strength  Curves  for  at  high,  nominal,  and  low  values,  with 
all  other  parameters  nominal,  and  <^  =  45  degrees. 


Figure  3.10:  Scattering  Strength  Curves  for  Sy  at  high,  nominal,  and  low  values,  with 
all  other  parameters  nominal,  and  (^  =  45  degrees. 


Figure  3.11:  Scattering  Strength  Curves  for  both  and  6y  at  high,  nominal,  and  low 
values,  with  all  other  parameters  nominal,  and  =  0  degrees. 


Chapter  4 

Model  Parameter  Estimation  from 
Monostatic  Data 


4.1  Parameter  Estimation  Methods 

As  discussed  in  Chapter  3,  it  is  desired  to  estimate  the  parameters  /i,  nm,  and 
S  given  backscatter  strength,  rntj,  incident  angle,  0,  and  acoustic  wavenumber,  k. 
The  data  that  will  be  used  for  the  estimation  is  mi,,  (0),  backscatter  strength  as  a 
function  of  incident  angle,  and  k  will  be  assumed  constant  for  this  vector.  The  vector 
^b3  (0)  will  typically  be  obtained  from  the  backscatter  strength  and  bathymetric  data 
produced  by  the  sonar  systems  mentioned  earlier.  A  single  vector,  and  correspond¬ 
ing  wavenumber,  is  obtained  for  each  transmit/receive  cycle  of  these  sonar  systems. 
The  work  discussed  in  this  paper  demonstrates  the  ability  to  estimate  the  desired 
parameters  by  generating  artificial  backscatter  data  using  the  BISSM  equation. 

The  estimation  of  //,  nm,  and  S  can  be  viewed  as  a  minimization  problem.  We  have 
mi„(0i),  the  scattering  strength  at  angle  0,  for  f  =  1  to  n.  We  also  have  a  function  T, 
{T  :  T(0i,j)},  which  is  the  BISSM  equation  and  is  assumed  to  be  a  model  for  nib,. 
In  T,  7  is  a  vector  composed  of  the  model  parameters:  7  =  [m,/z,^].  Thus  given  0, 
and  an  estimate  of  7  we  can  compute  an  estimate  for  scattering  strength,  mbs{0i)e- 
The  mean  square  error  between  the  estimate  and  the  actual  scattering  strength  is  a 
function  of  7,  and  is  denoted  £'(7).  We  want  to  minimize  this  error,  which  is  given 
by: 

£(7)  =  -D™*,(«.)-m.7)P  (4.1) 

".=1 

where  n  is  the  total  number  of  incident  angles.  By  minimizing  £(7),  our  estimate  of 
7  approaches  the  actual  value  70. 

Thus,  given  initial  estimates  of  these  parameters,  the  BISSM  equation  is  used  to 
generate  an  estimate  of  the  backscatter  strength  m(„(0)^.  A  better  estimate  of  the 
parameters  is  obtained  by  iteratively  adjusting  them  to  minimize  the  mean  square 
error  between  the  actual  backscatter,  mb,  {0),  and  the  estimated  backscatter  mj,  (0)^. 
Parameter  adjustment  is  typically  accomplished  through  random  search  or  gradient 
techniques.  While  random  search  techniques  can  guarantee  a  global  minimum  error. 
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they  are  very  computationaly  expensive.  There  are  several  gradient  search  methods 
available  for  nonlinear,  multidimensional  minimization  problems.  Techniques  that 
don’t  require  computation  of  first  derivatives  include  the  downhill  simplex  method, 
due  to  Nelder  and  Mead  [2],  and  Powell’s  method  [3].  Methods  that  require  the  com¬ 
putation  of  first  derivatives  include  the  Polak-Ribiere  algorithm  [4]  and  the  Broyden- 
Fletcher-Goldfarb-Shanno  algorithm  [5].  Gradient  techniques  can  provide  a  global 
minimum  if  the  error  surface  is  relatively  smooth.  For  highly  irregular  error  surfaces, 
which  are  encountered  with  noisy  data,  it  may  be  required  to  start  the  estimation  with 
multiple  initial  conditions  to  ensure  that  the  final  solution  is  not  a  local  minimum. 

Gradient  approaches  are  also  fairly  computationally  intensive,  and  a  gradient 
search  must  be  performed  for  each  new  data  vector  Tnbs(O)-  Another  approach  for 
the  estimation  of  7  is  to  use  a  feedforward  neural  network  to  model  the  inverse  rela¬ 
tionship,  F,  between  the  model  T,  and  7  where 

^(m7))  =  7  (4.2) 

Assuming  that  T  is  a  good  model  for  mbs  then  F{mt,s)  will  provide  a  good  estimate 
of  7.  As  shown  by  Kolmorogov  [6],  a  2  hidden  layer  feedforward  neural  network  can 
approximate  any  nonlinear  to  mapping  function  arbitrarily  close,  depending 
on  the  number  of  nodes  in  the  hidden  layers.  Furthermore,  neural  networks  are 
robust  in  the  presence  of  noise.  A  significant  advantage  of  a  neural  network  approach 
over  gradient  search  is  that  once  the  network  is  trained,  it  will  directly  produce 
an  estimate  of  7,  without  having  to  perform  a  gradient  search  for  each  new  mbs 
vector.  Also,  neural  networks  are  well  suited  for  implementation  on  parallel  machines. 
The  following  section  provides  a  brief  review  of  feedforward  neural  networks  and  the 
backpropagation  algorithm,  which  will  be  used  for  the  estimation  of  7  in  this  work. 


4.2  Review  of  Neural  Networks 

A  neural  network  is  a  device  that  can  be  used  to  recognize  signal  phenomena  or 
perform  numeric  functions.  A  neural  network  is  composed  of  one  or  more  layers,  each 
containing  one  or  more  neurons.  Each  neuron  in  a  feedforward  type  network  typically 
accepts  multiple  inputs,  where  these  inputs  are  signals  provided  to  the  network  or 
are  the  outputs  of  neurons  in  a  previous  layer.  A  single  neuron  is  essentially  a  linear 
combiner;  it  produces  a  weighted  sum  of  it  inputs.  Additionally,  each  neuron's  output 
is  limited  by  a  signum  or  sigmoid  function.  Signum  functions,  or  hard  limiters, 
are  typically  used  for  applications  where  the  network  is  intended  to  make  discrete 
classification  decisions  b2ised  on  the  inputs.  Sigmoid  functions,  or  soft  limiters,  are 
used  for  networks  designed  to  perform  various  real  valued  signal  processing  tasks. 

A  neural  network  typically  must  be  trained  to  perform  a  particular  function. 
Training  involves  providing  the  network  with  a  series  of  input  and  desired  signals, 
and  adapting  the  weights  of  the  neurons  within  the  network  to  minimize  the  error 
between  the  network’s  output  and  the  desired  signal.  One  of  the  more  popular  training 
algorithms  is  backpropagation  [7],  and  this  training  method  is  reviewed  here. 


1.5 


Updating  network  weights  with  the  backpropagation  algorithm  is  analogous  to 
the  process  used  by  the  adaptive  least  mean  square  (LMS)  filter.  The  algorithm 
computes  an  estimate  of  the  gradient  of  the  mean  square  error  with  respect  to  the 
system’s  weights,  and  this  information  is  used  to  move  the  weights  so  that  the  error 
approaches  a  minimum.  It  is  important  to  note  that  a  neural  network  is  a  nonlinear 
device,  typically  with  many  local  minima.  The  existence  of  local  minima  often  make 
it  difficult  for  the  network  to  reach  its  global  minimum  with  unsupervised  training. 
A  derivation  of  the  backpropagation  algorithm  is  now  presented. 

Given  Xk  as  the  input  vector  at  time  k,  and  w*.  as  the  weight  vector  at  time  k, 
the  output  of  the  linear  summation  of  a  single  neuron  is  given  by: 

Sk  =  xlv/k  (4. .3) 

and  the  output  of  the  neuron  is: 

yk  =  sgm(si)  (4.4) 


where  sgm  denotes  a  sigmoid  function.  The  hyperbolic  tangent  is  the  sigmoid  function 
used  in  this  paper.  The  error  here  is  defined  as  ejt  =  djt  —  yk  where  dk  is  the  desired 
signal  at  time  k. 

The  gradient  of  the  mean  square  error  with  respect  to  the  neuron's  weights  is 
given  by: 

where  E  denotes  the  expected  value.  By  eliminating  the  expected  value  we  obtain  a 
stochastic  estimate  of  the  gradient: 

|d-  (4.6) 

dwk 

resulting  in  the  following  weight  update  equation: 


WA..+  1  =  Wfc  -  j.1- 


del 


dwk 


(4.7) 


where  y.  is  the  learning  gain  that  controls  the  speed  of  convergence.  Taking  the  partial 
derivative  we  obtain  the  following: 

dck 


del 

dwk 


2ek 


2ek 


dwk 

d{dk  -  sgm(s^.)) 


(4.S) 


dsk 


dwk 

-2e*:Sgm'(sfc)- 

dwk 

-2ekSgm{sk)Xk 


However,  the  sigmoid  function  used  here  is  tanh,  so  we  have 


sgm'(.Sfc) 


d  tanh(s;t) 

d  Sk 

1  —  tanh^(s/t) 

I  -  vl 


(4.9) 
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(4.10) 


and  thus  the  weight  update  equation  for  a  single  neuron  becomes: 

wt+i  =  Wfc  +  Ofick  (l  -  yl)  Xk 

In  a  multilayer  neural  network,  desired  signals  are  typically  available  only  for  the 
neurons  in  the  output  layers.  A  general  weight  update  equation  for  any  neuron  is 
given  by  Widrow  and  Lehr  [7]  as: 


Wfc+i  =  Wfc  +  2fiSkXk  (4.11) 

where  6k  is  called  the  square  error  derivative  for  the  particular  neuron,  and  Xk  is  the 
input  vector  for  that  neuron.  For  the  output  layer  neurons,  where  desired  signals  are 
available,  6k  in  (4.11)  is  obtained  from  (4.10)  and  is  given  by  6k  =  6^.(1  —  yl).  For  a 
neuron  not  in  the  output  layer  6k  is  given  by  [7]: 

^L,k  =  (4-12) 

where  is  the  square  error  derivative  at  time  k  for  the  neuron  in  the  layer, 
s‘^  i^  is  the  output  of  the  same  neuron  (before  the  sgm  function)  at  time  k,  and 
is  the  backpropagated  error  for  that  neuron  at  time  k.  The  backpropagated  error  is 
given  by: 

^‘m.k  =  Z  ^IV^itn.n-^r.k  (4-13) 

»=1 

where  ^  is  the  weight  that  connects  the  neuron  in  the  layer  to  the 

neuron  in  the  (/  +  1)^‘  layer.  Thus  the  backpropagated  error  for  a  neuron  that  is  not 
in  the  output  layer  is  given  by  the  sum  of  the  square  error  derivatives  of  the  neurons 
in  the  following  layer,  each  scaled  by  the  neuron  weight  that  connects  the  neuron 
being  evaluated  with  the  neuron  in  the  next  layer. 
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Chapter  5 
Results 


In  this  cliapter  the  results  for  estimation  of  each  of  the  three  parameters  in  7 
are  shown,  using  noise-free  artificially  generated  data.  To  perform  the  estimation 
a  neural  network  program  was  written  in  ‘C’  using  the  backpropagation  algorithm 
[7].  This  program  is  given  in  Appendix  A,  and  the  program  that  controls  it  during 
training  is  in  Appendix  B.  The  program  implements  a  two  hidden  layer  network  with 
a  single  output  neuron,  and  provides  the  capability  to  select  the  number  of  neurons 
in  the  hidden  layers,  and  the  number  of  inputs  to  the  first  hidden  layer.  The  number 
of  neurons  in  each  of  the  two  hidden  layers  is  identical,  and  full  connectivity  between 
layers  has  been  used  thus  f  ir.  A  weight  jogging  capability  has  also  been  included, 
which  allows  the  addition  of  small  amounts  of  random  noise  to  the  weights  to  'jog’ 
the  network  out  of  local  minima.  This  feature  is  paramount  since,  as  illustrated 
by  Widrow  and  Lehr  [7],  the  error  surface  of  a  neural  network  is  typically  rich  in 
local  minima.  The  computation  time  has  been  found  to  be  a  nearly  linear  function 
of  the  total  number  of  weights  in  the  network.  With  learning  on  the  computation 
time  is  about  4.56  •  10“'^  seconds/iteration-tap  on  an  AT  computer.  With  learning  off 
the  time  is  about  one  third  of  this  value,  thus  the  backpropagation  process  requires 
about  two  thirds  of  the  total  computation  time.  Tests  on  a  Sun  4  computer  yielded  a 
computation  time  of  2.13  •  10"^  seconds/iteration  tap,  approximately  21  times  as  fast 
as  the  AT. 

The  network  architecture  that  is  best  for  a  given  problem  must  be  determined 
empirically,  typically  through  iterative  training  and  testing  of  different  architectures. 
The  method  used  here  is  to  start  with  a  small  network,  and  to  increase  its  size  to 
the  point  where  a  significant  reduction  in  error  is  no  longer  obtained.  The  data 
sets  for  training  are  generated  by  randomly  varying  the  parameters  in  the  BISSM 
equation,  providing  very  large  training  sets.  The  BISSM  equation  is  implemented  in 
the  program  given  in  Appendix  C.  The  training  method  used  here  starts  with  a  large 
learning  gain  for  faist  convergence,  followed  by  successive  reductions  in  gain  to  achieve 
lower  error.  Further  estimation  improvement  may  be  possible  by  adding  more  layers, 
or  through  advanced  techniques  such  as  momentum,  feedback,  or  feedforward  from 
nonadjacent  layers. 

The  following  sections  discuss  the  architecture,  training,  and  results  for  estimation 
of  p,  m,  and  S.  With  the  tanh  function  a  neuron’s  output  is  limited  to  the  range  of 
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[—1,1],  so  the  data  provided  to  the  network  is  scaled  and  in  some  cases  normalized 
to  improve  the  network’s  estimation.  The  range  of  values  for  p  and  6  are  sufficiently 
small  and  do  not  require  scaling.  The  values  for  n  and  m  may  be  large,  so  the 
network  is  trained  to  estimate  one-fourth  of  their  value.  As  previously  discussed, 
a  linear  relationship  is  a.ssumed  between  n  and  m.  The  relationship  used  here  is 
n  =  m  *  .3833  +  .51,  allowing  for  the  full  range  of  both  parameters.  Initial  training 
attempts  held  two  of  the  parameters  constant  and  randomly  varied  the  parameter  of 
interest.  It  was  found  that  the  error  could  be  made  small  for  the  parameter  being 
estimated,  but  could  become  large  if  the  other  two  parameters  were  varied  from 
the  values  used  during  training.  This  is  a  consequence  of  some  correlation  existing 
between  the  three  parameters.  To  reduce  this  problem,  all  three  parameters  are  varied 
randomly  throughout  their  full  ranges  during  training. 

5.1  Mackenzie  Coefficient  (fi)  Estimation 


Table  5.1:  p  network  training 


Learning 

Gain 

Iterations 

.Average  Percent 
Error 

0.8 

500 

136.4 

0.8 

500 

84.6 

0.6 

600 

42.5 

0.6 

600 

20.3 

0.4 

800 

7.75 

0.4 

800 

4.16 

0.2 

1000 

1.89 

0.2 

1000 

1.30 

0.1 

800 

0.87 

0.1 

800 

0.72 

0.05 

800 

0.58 

0.05 

800 

0.52 

0.01 

800 

0.50 

0 

800 

0.47 

For  the  estimation  of  the  input  vector  given  to  the  network  is  the  backscatter 
strength  for  incident  angles  ranging  from  15  to  60  degrees  in  2  degree  increments, 
where  each  input  to  the  network  corresponds  to  a  specific  incident  angle.  Angles 
higher  than  60  degrees  are  not  required  since  the  dominant  effect  of  /i  is  at  lower 
angles.  Angles  below  15  degrees  cannot  be  used  since  the  incident  angle  must  be 
large  enough  with  respect  to  n  in  the  Rayleigh  reflection  coefficient  calculation  to 
avoid  a  complex  result.  An  angle  increment  of  2  degrees  was  found  to  be  sufficient 
to  obtain  small  estimation  errors.  The  error  could  possibly  be  reduced  by  using  a 
smaller  increment,  but  this  will  also  require  a  larger  network.  Since  the  backscatter 
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Figure  5.1;  Probability  density  function  of  the  percent  estimation  error  for  the 
estimating  program.  The  scattering  strength  data  was  obtained  by  randomly  varying 
S,  /i,  and  m. 

strength  is  small  at  lower  angles  (ranging  from  about  10~®  to  10“^),  the  data  vector  is 
multiplied  by  10  before  giving  it  to  the  network.  Also,  the  data  vector  was  ‘centered’ 
by  subtracting  0.15  from  all  values.  This  allows  the  data  at  the  lower  angles  to  have 
approximately  equal  influence  on  the  training  process  as  the  higher  angles. 

The  final  configuration  for  the  n  network  has  23  inputs  to  each  of  the  nodes  in 
the  first  layer  and  10  nodes  in  each  of  the  two  hidden  layers.  The  n  estimation 
program  is  given  in  Appendix  D.l.  The  final  training  sequence  is  shown  in  Table 
5.1,  and  an  average  percent  error  of  0.47  Wtis  obtained  with  training  off  (gain  =  0). 
As  shown  in  this  table  the  network  converged  very  rapidly  to  a  small  error.  Faster 
computation  rates  can  be  obtained  by  increatsing  the  angle  increment  in  the  input  data 
vector,  probably  with  a  minimal  increase  in  estimation  error.  The  probability  density 
function  of  the  percent  estimation  error  is  shown  in  Figure  5.1,  and  indicates  that  for 
the  majority  of  the  test  data  the  estimation  error  is  below  1  percent.  As  expected,  the 
estimation  of  fi  was  relatively  easy  since  its  dominant  effect  on  backscatter  strength 
is  essentially  a  near  uniform  change  in  amplitude  at  all  lower  incident  angles. 


5.2  m  Estimation 

For  the  estimation  of  m  the  input  vector  given  to  the  network  is  the  backscatter 
strength  for  incident  angles  ranging  from  71  to  88  degrees  in  1  degree  increments, 
where  each  input  to  the  network  corresponds  to  a  specific  incident  angle.  Angles 
lower  than  71  degrees  are  not  required  since  the  dominant  effect  of  nm  is  at  higher 
angles,  and  exclusion  of  the  lower  angles  helps  to  reduce  sensitivity  of  the  estimation 
to  changes  in  /x.  Angles  of  89  and  90  degrees  gave  backscatter  strength  values  greater 
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Table  5.2:  m  network  training 


Learning 

Gain 

Iterations 

Average  Percent 
Error 

1 

2000 

4.55 

0.9 

2000 

1.78 

0.8 

2000 

1.45 

0.7 

2000 

1.35 

0.6 

2000 

1.19 

0.5 

2000 

1.10 

0.4 

2000 

1.03 

0.3 

2000 

1.02 

0.2 

2000 

0.98 

0.1 

2000 

0.91 

0.05 

2000 

0.88 

0.0 

1000 

0.87 

than  1  and  were  excluded.  An  angle  increment  of  1  degree  was  found  to  be  sufficient 
to  obtain  small  estimation  errors.  In  the  m  network  the  parameter  being  estimated  is 
m,  but  as  mentioned  previously  a  linear  relationship  is  being  used  to  obtain  n  given 
a  value  for  m.  The  error  could  possibly  be  reduced  by  using  a  smaller  increment, 
but  this  will  also  require  a  larger  network.  The  data  vector  is  normalized  to  reduce 
the  effect  of  changes  in  fi  by  subtracting  the  backscatter  strength  at  the  angle  of  71 
degrees  from  all  other  points  in  the  vector.  Also,  the  data  point  corresponding  to  71 
degrees  in  the  data  vector  is  set  to  1  to  provide  the  network  with  an  adjustable  offset 
capability.  As  previously  mentioned,  the  values  of  m  may  become  too  large  for  the 
tanh  function,  so  the  network  is  trained  to  estimate  m/4.  Thus  the  final  estimate  of 
m  is  given  by  the  network  output  multiplied  by  4. 

The  final  configuration  for  the  m  network  has  18  inputs  to  each  of  the  nodes 
in  the  first  layer  and  8  nodes  in  each  of  the  two  hidden  layers.  The  m  estimation 
program  is  given  in  Appendix  D.2.  The  final  training  sequence  is  shown  in  Table 
5.2,  and  an  average  percent  error  of  0.87  was  obtained  with  training  off.  The  table 
indicates  that  the  network  trained  rapidly  to  an  error  of  about  5  percent,  but  then 
slowly  decreased  to  its  final  value.  This  indicates  that  the  error  surface  is  shallow  in 
the  vicinity  of  the  minimum,  which  will  cause  difficulty  in  training  with  noisy  signals. 
Recall  that  all  three  parameters  (/i,  m,  and  6)  are  being  varied  randomly,  so  this 
network  has  successfully  differentiated  between  the  effect  of  m  and  S  on  the  slope  of 
the  backscatter  strength  curve  at  high  angles.  The  probability  density  function  of  the 
percent  estimation  error  is  shown  in  Figure  5.2,  and  indicates  that  for  the  majority 
of  the  test  data  the  estimation  error  is  below  2  percent. 
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Figure  5.2;  Probability  density  function  of  the  percent  estimation  error  for  the  m 
estimating  program.  The  scattering  strength  data  was  obtained  by  randomly  varying 
S,  /i,  and  m. 

5.3  6  Estimation 

For  the  estimation  of  6  the  input  vector  given  to  the  network  is  the  backscatter 
strength  for  incident  angles  ranging  from  71  to  88  degrees,  and  an  increment  of  0.5 
degrees  was  required  to  obtain  low  estimation  errors.  As  with  the  m  parameter, 
angles  lower  than  71  degrees  are  not  required  since  the  dominant  effect  of  ('’  is  at 
higher  angles,  and  exclusion  of  the  lower  angles  helps  to  reduce  sensitivity  of  the 
estimation  to  changes  in  fi.  The  data  vector  for  the  6  network  is  also  normalized 
using  the  71  degree  data  point,  and  this  point  is  set  to  1  to  provide  the  network  with 
an  adjustable  offset  capability. 

The  final  configuration  for  the  S  network  has  36  inputs  to  each  of  the  nodes  in  the 
first  layer,  and  7  nodes  in  each  of  the  two  hidden  layers.  The  S  estimation  program 
is  given  in  Appendix  D.3.  The  final  training  sequence  is  shown  in  Table  5.3,  and 
an  average  percent  error  of  5.16  was  obtained  with  training  off.  Training  of  the  6 
network  proved  to  be  more  difficult  than  for  m  and  /i.  The  training  sequence  shows 
rapid  convergence  to  an  error  of  about  20  percent,  but  a  much  slower  convergence  to 
smaller  errors  indicating  a  very  flat  error  surface  near  the  minimum.  Using  a  network 
with  only  IS  inputs  to  the  first  layer  yielded  a  minimum  error  of  about  15  percent 
without  the  adjustable  offset  capability  (71  degree  data  point  =  1),  and  about  12 
percent  with  this  capability.  In  earlier  tests,  when  /y  and  m  were  held  constant,  an 
error  of  1  percent  was  achieved,  but  this  became  as  large  as  20  percent  when  different 
values  of  /r  and  nm  were  used  in  testing  than  those  used  during  training.  Figure 
5.3  shows  the  probability  density  function  of  the  percent  estimation  error  for  the  S 
network,  and  indicates  that  for  the  majority  of  the  test  data  the  estimation  error  is 
below  about  10  percent. 
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Table  5.3:  8  network  training 


Learning 

Gain 

Iterations 

Average  Percent 
Error 

1.6 

4000 

17.1 

1.5 

4000 

12.6 

1.4 

4000 

8.74 

1.3 

4000 

7.89 

1.2 

4000 

7.70 

1.1 

4000 

7.10 

1.0 

4000 

6.98 

0.9 

4000 

7.05 

0.8 

4000 

6.72 

0.7 

4000 

6.50 

0.6 

4000 

6.07 

0.5 

4000 

6.03 

0.4 

4000 

6.09 

0.3 

4000 

5.99 

0.2 

4000 

5.80 

0.1 

4000 

5.53 

0.05 

4000 

5.51 

0.0 

1000 

5  16 

Pcrceni  Error 


Figure  5.3:  Probability  density  function  of  the  percent  estimation  error  for  the  8 
estimating  program.  The  scattering  strength  data  wa.s  obtained  by  randomly  varying 
(5,  /i,  and  m. 


Chapter  6 
Summary 


This  work  has  demonstrated  the  ability  to  estimate  certain  parameters  in  the 
BISSM  [1]  acoustic  scattering  model  using  artificial  backscatter  data.  The  six  param¬ 
eters  of  interest  are  /x,  n,  m,  tr,  61,  and  8y.  It  was  desired  to  estimate  these  parameters 
given  inputs  of  backscatter  strength,  nibs,  the  incident  angle,  0,  and  the  azimuth  angle, 
(j>.  These  inputs  potentially  can  be  obtained  from  backscatter  and  bathymetric  data 
collected  by  sidescan  or  multibeam  acoustic  imaging  systems.  Of  the  six  parameters 
only  three  could  be  estimated;  n  and  m  are  highly  correlated,  cr  has  minimal  effect 
on  the  backscatter  strength  at  1.2  kHz,  and  Sx  and  8y  are  highly  correlated.  Sx  and 
were  assumed  to  be  equal,  giving  a  single  parameter  S,  and  this  eliminated  6  from 
the  monostatic  version  of  the  BISSM  equation.  A  linear  relationship  was  assumed 
between  n  and  m,  allowing  m  or  n  to  be  independently  estimated.  Feedforward  neural 
networks  were  used  to  estimate  the  three  parameters,  /x,  m,  and  6,  given  backscatter 
strength  values  as  a  function  of  the  incident  angle.  Using  backscatter  data  generated 
with  the  BISSM  equation,  the  networks  successfully  estimated  /x  and  in  with  less  than 
1  percent  average  error,  and  8  with  about  5  percent  average  error. 

Further  developmen.  of  the  neural  networks  used  in  this  project  will  be  required  for 
application  to  real  data.  Real  data  will  necessarily  be  corrupted  with  noise,  and  will 
seldom  have  ground  truth  information  available.  Without  ground  truth,  real  data 
cannot  be  used  to  train  the  networks.  One  approach  to  improving  the  parameter 
estimates  in  the  presence  of  noise  would  be  to  analyze  the  noise  character  in  real 
data.  Then,  synthetic  noise  of  similar  character  can  be  used  to  contaminate  artificially 
generated  data,  and  the  networks  can  be  retrained  to  improve  their  performance  with 
noisy  data.  Also,  the  number  and  range  of  incident  angles  available  from  real  data 
depends  on  the  bottom  morphology  as  well  as  the  survey  system,  with  the  result  that 
the  networks  may  not  have  a  full  input  data  vector.  Further  testing  is  required  to 
determine  the  effect  on  the  network’s  estimates  in  the  event  of  missing  input  data 
points. 
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Appendix  A 

Neural  Network  Subroutine 


/♦  nn .  c 

Brian  Bourgeois  Created:  26AUG91  Last  Mod:  26AUG91 
This  neural  net  sub  is  designed  to  take  a  fixed  vector 
length  data  as  input  and  produce  a  single  output  variable, 
representing  some  mapping  from  R“n  to  R“l.  Full  connectivity 
is  used,  and  backpropagation  is  the  learning 
algorithm.  The  ability  to  ‘jar’  the  weights  during  learning 
is  included  (annealing) .  A  call  to  nn  causes  a  single  pass 
through  the  network. 

All  controlling  information  is  passed  through  vectors  arch 
(architecture)  and  train.  Data  is  passed  thru  in,  des,  and  est . 
in  is  the  input  vector,  des  is  the  desired  vector,  and  est  is 
the  estimate  vector.  The  network  weights  are  passed  via  taps. 

In  the  vector  arch,  it  is  possible  to  specify  the  number 
on  inputs  to  each  of  the  first  layer  nodes,  and  the  number  of 
nodes  in  each  layer.  A  variable  for  specifying  the  number  of 
layers  is  used,  but  the  program  is  only  set  up  for  3  layers 
at  this  time. 

In  the  vector  train,  the  learning  (adaptive)  gain  is 
specified,  which  is  typically  less  than  1.  The  learning  switch 
controls  whether  the  node  errors  are  computed  and  weights  updated. 

The  annealing  gain  controls  the  eonount  of  noise  added  to  the 
weights  in  a  iteration,  train [3]  and  train [4]  are  not  used  in 
this  sub,  but  in  the  controlling  program. 

*/ 


/*  train [0]  learning  gain  */ 
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/*  train [l]  learning  switch,  1  =  on  */ 

/*  train [2]  annealing  gain,  0  =  off  ♦/ 

/♦  train[3]  0  =  init  taps,  1  =  load  from  taps. in  */ 

/*  train [4]  Number  of  repetitions  for  this  set  */ 

/♦  arch[0]  =  Number  of  layers  */ 

/*  arch[l]  =  Number  of  neurons  in  first  layer  */ 

/♦  arch [2]  =  Number  of  neurons  in  second  layer  */ 

/*  arch [3]  =  Number  of  neurons  in  output  layer  */ 

/*  arch[4]  =  Number  of  inputs  to  first  layer  */ 

#include  "nn.h" 

/4>«******4>***«*4<**«****«******«***4c«*:«'****1c*«***>t'!ti*/ 

int  nn(in,des,est , arch, train, taps) 

double  *in;  /*  in  =  input  vector  ♦/ 

double  *des;  /*  des  =  desired  vector  *! 

double  ♦est ;  /♦  est  =  estimated  vector  */ 

long  *arch;  /*  network  configuration  parajneters  */ 

double  *train;  /*  network  training  paraoneters  */ 

double  ♦**taps;  /*  network  weights  */ 

{ 

/♦♦♦♦**  declare  variables  ♦♦♦♦♦*♦/ 

/*  network  configuration  ♦/ 

/♦  arch[0]  =  number  of  layers  of  neurons  ♦/ 

/*  arch[l]  =  layer  1  length  */ 

/*  arch [2]  =  layer  2  length  */ 

/♦  arch [3]  =  output  layer  length  */ 

/*  arch [4]  =  Number  of  inputs  to  first  layer  */ 

/*  train [0]  learning  gain  */ 

/♦  train [l]  learning  switch,  1  =  on  ♦/ 

/♦  train [2]  annealing  gain,  0  =  off  */ 

/*  train [3]  0  *  init  taps,  1  =  load  from  taps. in  */ 

/♦  train [4]  Number  of  repetitions  for  this  set  */ 

long  no_layers,  /*  Number  of  network  layers  with  neurons  *! 
*no_neurons ,/♦  Number  of  neurons  in  each  layer  */ 
*no_inputs,  /♦  Number  of  inputs  for  neuron  in  given  layer*/ 
layer,  /*  Current  layer  number  */ 

node,  /*  Current  node  number  */ 

tap,  /*  Current  tap  number  in  a  neuron  */ 
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double 


elem; 

sum, 

**node_out , 

**err, 

corr; 


/*  Node  in  next  layer,  for  error  calc  */ 
/♦  Summer  output  for  a  neuron  */ 

/♦  Outputs  on  each  node  in  nn  */ 

/*  Error  for  each  neuron  */ 

/*  Tap  correction  for  updata  */ 


/♦  Misc  variables  */ 
long  Ctrl,  /*  counter  */ 
ctr2,  /*  counter  ♦/ 
ctr3;  /*  counter  */ 
double  temp: 


/♦♦**♦*  Load  up  architecture  variables  ******/ 
no.layers  =  arch[0]; 


no.neurons  =  (long  *)malloc(sizeof (long)*(int)arch[0] ) ; 
if (no_neurons  ==  NULL){ 

printf( "no .neurons  allocation  error  {nn}\n"); 
return: 

} 


no.neurons [0]  =  arch  Cl] 
no_neurons[l]  =  arch [2] 
no_neuronsC2]  =  arch [3] 


no. inputs  =  (long  *)malloc(si2eof (long)*(int)arch[0] ) ; 
if (no. inputs  ==  NULL){ 

printf ("no. inputs  allocation  error  {nn}\n") ; 
return: 

} 

no.inputsCO]  =  arch[4]:  /*  Inputs  to  first  layer  */ 

no. inputs [l]  =  no. neurons [0] :  /♦  Inputs  to  second  layer  */ 

/*  same  as  no. neurons  in  first  layer  for  full  connectivity  */ 
no_inputs[2]  =  no.neuronsCl] :  /*  Inputs  to  output  layer  */ 

/*  same  as  no.neurons  in  second  layer  for  full  connectivity  */ 

/******  Memory  allocation  ♦***♦♦/ 

node. out  *  (double  ♦♦)malloc(sizeof (double  **) *no. layers) ; 
if (node. out  ==  NULL){ 

printf ("node. out  allocation  error,  level  l\n") ; 
return; 

}  /*  if  node. out  *=  NULL  */ 


f or (ctr 1=0; ctrl<no. layers ;ctrl++){ 

node. out [ctrl]  =  (double  *)malloc(sizeof (double) *no_neurons[ctrl] ) 
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if (node.out [ctrl]  ==  NULL){ 

printf ("node_out  allocation  error,  level  2\n") : 
return; 

}  /♦  if  node.out []  ==  NULL  ♦/ 

}  /♦  for  Ctrl  ♦/ 

err  *  (double  **)malloc(sizeof (double  ♦*)*no_layers) ; 
if (err  «  NULL){ 

printf ("err  allocation  error,  level  l\n") ; 
return ; 

}  /♦  if  err  ==  NULL  */ 
for(ctrl=0; ctrl<no_layers:ctrl++){ 

errCctrl]  =  (double  *)malloc(sizeof (double)*no_neurons [ctrl] ) ; 
if(err[ctrl]  ==  NULL){ 

printf ("err  allocation  error,  level  2\n") ; 
return; 

}  /♦  if  err[]  ==  NULL  */ 

}  /*  for  Ctrl  */ 


Z********:^***  Functional  Part  of  Program  *♦♦♦♦***♦**♦*/ 

Forward  Sweep  through  network  *****/ 

/*  Layer  1  */ 
layer  *  0; 

f  or (node*0 ; node<no_neurons [layer] ; node++) { 
sum*0 ; 

f or (tap=0 ;tap<no_ inputs [layer] ;tap++) { 

sura  +=  taps [layer] [node] [tap]* in [tap] ; 

}  /*  for  tap  */ 

node_out  [layer]  [node]  =  tanh(suin) ; 

}  /*  for  node  */ 

/♦  Following  Layers  */ 
f or (layer*l ; layer<no_layers ; layer++) { 

f or(node=0 ;node<no_neurons [layer] ;node++){ 
sum=0 ; 

f or(tap=0;tap<no_inputs [layer] ;tap++){ 

sum  +=  taps [layer] [node] [tap] *node_out [layer- l] [tap] ; 
}  /*  for  tap  */ 

node_out  [layer]  [node]  =  tanh(suin)  ; 

}  /*  for  node  ♦/ 

}  /♦  for  layer  */ 
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/♦  Estimated  Outputs  */ 

f or (ctrl=0 ; ctr l<arch [3] ; ctrl++) { 

estCctrl]  =  node.out [no_layers-l] [ctrl] ; 

} 

/♦  printf  ("train [1]  =  V.fXn"  .train [l]  )  ;  */ 

if(train[l]  ==  !){  /*  Do  this  next  section  for  learning  on  only  */ 

/♦♦**♦  Backward  Sweep  through  network  ***♦*/ 

/♦  Compute  Errors  */ 

/*  Output  Layer  */ 
layer  =  no.layers-1; 

f or (node=0 ; node<no_neurons [layer] ; node++) { 

err [layer] [node]  =  (des [node] -est [node] ) * ( 1-est [node] *est [node] ) ; 

} 

/♦  Following  Layers  */ 
f or(layer®no_layers-2 ; layer>-l ; layer — ) { 

for (nodeaO ; node<no_neurons [layer] ; node++) { 
sum  *  0; 

f or(elem=0 ;elem<no_neurons [layer+l] ; elem++){ 

sum  +=  err [layer+l] [elem]*taps [layer+l] [elem] [node] ; 

}  /*  for  elem  */ 

temp  =  node_out [layer] [node] *node_out [layer] [node] ; 
err [layer] [node]  =  sum* (1-terap) ; 

}  /♦  for  node  */ 

}  /*  for  layer  */ 

/♦  Update  filter  weights  */ 

/*  All  but  first  layer  */ 
f or (layer=no_layers-l ; layer>0 ; layer — ) { 

f or (node=0 ;node<no_neurons [layer] ; node++) { 
f or(tap=0 ;tap<no_inputs [layer] ;tap++){ 

corr  =  2.*train[0]*err[layer] [node]*node_out[layer-l] [tap] ; 
taps  [layer]  [node]  [tap]  +*  corr; 

}  /*  for  tap  */ 

>  /*  for  node  */ 

}  /♦  for  layer  */ 

/*  first  layer  */ 
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layer  =  0; 

f or(node=0 ;node<no_neurons [layer] ;node++) { 
f or (tap=0 ; tap<no_inputs [layer] ; tap++) { 

corr  =  2. *t rain [0] ♦err [layer]  [node] *111  [tap]  ; 
taps  [layer]  [node]  [tap]  +=  corrj 
}  /*  for  tap  */ 

}  /♦  for  node  */ 

/*  Jog  taps:  annealing  method  */ 
if (train[l]”l  k&  train[2]>0){ 

/*  printfC'Jog  taps\n");  */ 
f or (layer=0 ; layer<no_layers ; layer++) { 

f or (node=0 ;node<no_neurons [layer] ; node++) { 
for (tap=0;tap<no_inputs [layer] ;tap++){ 

taps [layer] [node] [tap]  +*  train[2]*(((double)rand()/RAND_MAX)- .5 
}  /*  for  tap  ♦/ 

}  /♦  for  node  */ 

}  /♦  for  layer  */ 

}  /♦  if  train [l]  &&  train [2]  ♦/ 

}  /♦  if  train [1]  ==  1,  Learning  on  ♦/ 

End  Sub  ♦♦♦♦♦♦/ 

/*  free  memory  ♦/ 

f or(ctrl=0; ctrl<no_layers ;ctrl++){ 
free(node_out [ctrl] ) ; 

} 

free(node_out) j 

f or(ctrl=0; ctrl<no_layers ;ctrl++){ 
free(err [ctrl] ) ; 

} 

free(err) ; 

f ree(no_neurons) ; 
free(no_inputs) ; 

return; 

> 
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Appendix  B 

Neural  Network  Training  Program 


/♦  del.c 

Brian  Bourgeois  Created:  13SEP91  Last  Mod:  23SEP91 
This  is  the  driving  main  program  for  the  neural  net  sub  'nn.c'. 
The  taps  array  is  created  and  maintained  in  this  program,  and 
network  architecture  and  training  motif  are  also  specified  here 
2Lnd  passed  to  the  network  sub.  Each  call  to  nn  does  a  single 
pass  through  the  network. 

This  prograun  also  provides  an  area  for  calling  a  data  generation 
routine,  for  specifying  the  input  and  desired  signals.  */ 

#include  "nn.h" 

♦include  "bism.h" 

extern  int  nn(double*,  double*,  double  *,  long  *, 
double  *,  double  ***) ; 

extern  int  bism(double  *ss,  double  *theta,  double  *model) ; 

/*  exclude  this  for  cc  */ 

mainCint  argc,  char  *argv[]) 

/*  main(argc,argv) 

int  argc; 

char  *argv  []  ;  ♦/ 

{ 

/*it^****  Declare  variables  ♦*♦*****/ 

/*  Data  vectors  */ 

double  *in;  /*  network  input  vector  */ 
double  *des;  /*  network  desired  vector  */ 


double  *est; 


/*  network  estimated  vector  */ 


/*  network  configuration  */ 
long  *arch;  /♦  network  configuration  parameters  */ 
double  ♦train;  /*  network  training  parajneters  ♦/ 
double  ♦♦♦taps;  /♦  network  weights  ♦/ 

long  no_layers;  /♦  Number  of  network  layers  with  neurons  ♦/ 

long  ♦no.neurons;  /♦  Number  of  neurons  in  each  layer  */ 

long  ♦no.inputs;  /♦  Number  of  inputs  for  a  neuron  in  a  given  layer  */ 

/♦  Training  variables  ♦/ 

long  num.sets;  /♦  Number  of  training  sets  in  train.dat  file  ♦/ 

long  set_cnt;  /♦  Current  training  set  ♦/ 

double  error;  /♦  Output  error  ♦/ 

double  mse;  /♦  Output  mse  ♦/ 

double  perr;  /♦  desired  signal  squared  ♦/ 

double  errdev;  /♦  error  std  deviation  ♦/ 

double  ♦iterr;  /♦  Error  for  each  iteration  ♦/ 

/♦  Misc  variables  ♦/ 


long  Ctrl,  /»  '.ounter  ♦/ 
ctr",  ^  counter  ♦/ 
c+rS,  /♦  counter  ♦/ 

^'lag. 

temp; 

FILL  ♦ftaps.in,  /♦  Input  taps  file,  with  arch  header  ♦/ 

♦ftaps.out,  /♦  Output  taps  file,  with  arch  header  */ 
♦ftrain,  /♦  File  with  training  instructions  ♦/ 
♦fdefarch,  /♦  default  network  architecture  file  ♦/ 
♦fdef model,  /♦  Default  model  paraimeters  file  ♦/ 

♦fiterr;  /♦  iteration  error  data  dump  file  ♦/ 

/♦  Model  variables  ♦/ 


double  ♦theta,  /♦  incident  angle  array  ♦/ 
♦ss,  /♦  scattering  strength  data  ♦/ 
model [9];  /♦  model  parameters  ♦/ 


/♦  Model 
model  [O] 
model [1] 
model [2] 
model [3] 
model [4] 


Parajneters 

=  size  of  data  vectors  theta  and  ss 
=  n.  Ratio  of  sound  speeds 

=  m.  Ratio  of  densities 

*  mu,  Lambert  coefficient 

=  sigma.  Microscale  heights  roughness 
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model [5]  *  phi,  azimuth  angle  in  radians 

model [6]  =  delx,  RMS  slope  deviation  along  track,  radians 

model [7]  =  dely,  RMS  slope  deviation  across  track,  radians 


model [8] 

=  k; 

Wavenumber 

*! 

long 

cntl , 

/* 

counter  */ 

cnt2 , 

/* 

counter  */ 

min. 

/♦ 

minimum  thet 

angle  */ 

max ; 

/* 

maximum  thet 

emgle  */ 

double  step; 

/* 

thet  angle  step  size  */ 

/*  Assign  model  parameters  */ 
flag  =  0; 

min  =71;  /*  minimum  theta  in  degrees  */ 

max  *  88;  /*  max  theta  :  60  to  89  for  del  test  */ 

step  =  .5;  /*  theta  angle  step  size  */ 

/*  compute  vector  length  */ 

modelCO]  =  (double) (max  -  min  +  l)/step  ;  /*  size  of  theta  vector  */ 
/*  Note;  No  of  Inputs  to  network  should  be  same  as  model [0]  */ 

Obtain  network  configuration  data  ******/ 

/♦  Allocate  storage  for  controlling  parameters  */ 

arch  =  (long  *)malloc(sizeof (long)*(int)5) ; 
if (arch  ==  NULL){ 

printf("arch  allocation  error\n"); 
return; 

} 

train  =  (double  *)malloc(sizeof (double) *(int)5) ; 
if (train  ==  NULL){ 
printf ("train  allocation  errorXn") ; 
return; 

} 

/♦  Pick  up  first  training  set.  This  is  done  first  to  determine 
(from  init  field)  if  the  architecture  will  be  specified  by  the  file 
defarch.dat,  or  if  weights  from  a  previous  execution  will  be  used, 
and  thus  the  arch  will  be  contained  in  a  file  taps. in  */ 

ftrain  =  fopen("train.dat","r") ; 
if (f train  ==  NULL){ 
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printf ("Training  file  does  not  exist  —  terminatingXn") ; 
return; 

} 

/♦  read  in  no  of  training  sets  ♦/ 
fsceuif  (ftrain,"  7, Id"  ,&nuin_sets)  ; 

/*  read  in  first  training  set  ♦/ 
f or(ctrl=0 ; ctrl<5 ;ctrl++){ 

fscanfCf train,"  */,lf  "  ,&train[ctrl3 )  ; 

}  /*  for  Ctrl  */ 

/♦  train [0]  learning  gain  */ 

/♦  train [1]  learning  switch,  1  =  on  */ 

/♦  train [2]  annealing  gain,  0  =  off  */ 

/♦  train [3]  0  =  init  taps,  1  =  load  from  taps. in  */ 

/*  train [4]  Number  of  repetitions  for  this  set  */ 

if(trainC3]  ==  l){ 

/*  Read  in  architecture  from  existing  taps  file,  taps. in  */ 
ftaps.in  =  f open("taps . in", "r") ; 
if(ftaps_in  ==  NULL){ 

printf ("taps . in  does  not  exist  —  terrainating\n") ; 
return; 

} 

for(ctrl=0; ctrl<5 ;ctrl++){ 

fscanf  (ftaps_in,"  ’/.Id"  ,&arch[ctrl]  ) ; 

}  /*  for  Ctrl  */ 

}  /♦  if  train [3]  ♦/ 
else{ 

/*  use  default  network  architecture  vector  ♦/ 
fdefarch  =  f open("defarch .dat" , "r") ; 
if(fdefarch  ==  NULL){ 

printf ("defarch.dat  does  not  exist  —  terminatingXn") ; 
return; 

> 

f or(ctrl®0; ctrl<5 ;ctrl++) { 

fscanf  (fdefarch,  "  '/.Id"  ,&arch[ctrl]  ) ; 

}  /♦  for  Ctrl  */ 
f close(f def arch) ; 

/*  arch[0]  =  Number  of  layers  */ 

/*  archCl]  *  Number  of  neurons  in  first  layer  */ 

/*  arch [2]  =  Number  of  neurons  in  second  layer  */ 

/*  arch [3]  =  Number  of  neurons  in  output  layer  */ 


/♦  arch [4]  *  Number  of  inputs  to  first  layer  */ 

}  /♦  else  */ 

/*  Load  up  architecture  variables  */ 
no_layers  =  arch[0] ; 

no_neurons  =  (long  *)malloc(sizeof (long)*(int)arch[0] ) ; 
if (no.neurons  ==  NULL){ 
printf ("no_neurons  allocation  error\n"); 
return; 

} 

no_neurons [O]  =  arch[l] ; 
no_neurons [l]  =  arch [2] ; 
no_neurons [2]  =  arch [3] ; 

no_inputs  =  (long  *)malloc(sizeof (long)*(int)arch[0] ) ; 
if(no_inputs  ==  NULL){ 
printf ("no_inputs  allocation  error\n"); 
return; 

> 

no_inputs[0]  =  arch [4] ;  /♦  Inputs  to  first  layer  ♦/ 
f*  Note  this  should  be  same  as  length  of  ss  vector  */ 
no_inputs[l]  =  no.neurons [0] ;  /♦  Inputs  to  second  layer  */ 

/*  same  as  no. neurons  in  first  layer  for  full  connectivity  */ 
no_inputs[2]  =  no.neurons [l] ;  /*  Inputs  to  output  layer  */ 

/*  same  as  no.neurons  in  second  layer  for  full  connectivity 

/*♦***♦  Allocate  storage  *♦♦4**/ 

in  =  (double  *)malloc(sizeof (double)*(int)arch [4] ) ; 
if (in  ==  NULL){ 

printf ("in  allocation  error\n") ; 
return; 

> 


des  =  (double  *)malloc(sizeof (double)*(int)arch[3] ) ; 
if(des  ==  NULL){ 

printf ("des  allocation  error\n"); 
return; 

} 

est  =  (double  *)malloc(sizeof (double)*(int)arch[3] ) ; 
if(est  ==  NULL){ 

printf ("est  allocation  error\n") ; 
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return: 

} 


theta  =  (double  *)malloc(sizeof (double)*(int)model [0] ) ; 
if  (theta  ==  NULL)-C 

printf ("theta  allocation  errorXn") ; 
return; 

} 


ss  *  (double  *)malloc(si2eof (double)*(int)model[0]  ) ; 
if(ss  ==  NULL)-C 

printf ("ss  allocation  error\n") ; 
return; 

} 

/*  nn  weights  allocation  is  a  bit  more  involved  */ 
taps  =  (double  ***)malloc(sizeof (double  ***) *no_layers) ; 
if (taps  ==  NULL){ 

printf ("taps  allocation  error,  level  l\n") ; 
return; 

}  /*  if  taps  ==  NULL  */ 
for(ctrl=0;ctrl<no_layers ;ctrl++){ 

taps [Ctrl]  *  (double  ♦*)malloc(sizeof (double  **)*no_neurons [ctrl] ) ; 
if (taps[ctrl]  ==  NULL){ 

printf ("taps  allocation  error,  level  2\n"); 
return; 

}  /*  if  taps[]  ==  NULL  */ 

}  /*  for  Ctrl  */ 

f or (ctrl*0 ; ctr l<no_layers ; ctrl++) { 
f or(ctr2=0; ctr2<no_neurons [ctrl] ;ctr2++){ 

taps [ctrl] [ctr2]  =  (double  *)malloc(sizeof (double) *no_inputs [ctrl] ) ; 
if  (taps  [ctrl]  [ctr2]  ==  NULLX 
printf ("taps  allocation  error,  level  3\n"); 
return; 

}  /*  if  taps[][]  ==  NULL  */ 

}  /*  for  ctr2  */ 

}  /*  for  Ctrl  */ 

/*  Working  Area  MAIN  */ 

/ Load  taps  array  ******/ 
if (train[3] ==!){ 


/*  Load  taps  from  taps. in  and  close  file  ♦/ 
f or(ctrl=0; ctrl<no_layers ;ctrl++){ 

f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 
f or (ctr3=0 ; ctr3<no_inputs [ctrl] ; ctr3++) { 
f  scanf  (ftaps_in, "  ‘/,le"  ,&taps  [ctrl]  [ctr2]  [ctr3]  ) ; 

}•  /*  for  tap  ♦/ 
y  /*  for  node  */ 

}  /♦  for  layer  */ 
f close(ftaps_in) ; 

}  /♦  if  train [3]  */ 
else-( 

/♦  Initialize  taps  with  small  random  numbers  */ 
f or (ctrl=0 ; ctrl<no_layers ; ctrl++) { 
f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 
f  or  (ctr3=0 ;  ctr3<no_ inputs  [ctrl]  ;ctr3++)  ■[ 

taps [ctrl] [ctr2] [ctr3]  =  .5*(((double)rand()/RAND_MAX)- .5) ; 

/*  -.5  to  .5  */ 

}  }  } 

/♦  Print  weights 

for(ctrl=0; ctrl<no_layers ;ctrl++){ 
f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 
f or (ctr3=0 ; ctr3<no_inputs [ctrl] ; ctr3++) { 
printf  ("'/if  \n",  taps  [ctrl]  [ctr2]  [ctr3])  ; 

}  }  }  */ 

}  /♦  else  */ 

/*  Assign  incident  amgles  */ 
f or(cntl=0; cntl<model [O] ; cntl++){ 

theta[cntl]  =  (double) (min  +  cntl*step)  *  PI  /  180.; 

/♦  printf  ("theta  =  y,f  \n" , theta [cntl]  )  ;  */ 

} 

/♦  Load  nominal  model  pareuneters  ♦/ 

fdefmodel  =  f openC'defmodel .dat" ,"r") ; 
if(fdefmodel  ==  NULL){ 

printf ("defmodel .dat  does  not  exist  —  terminatingXn") ; 
return; 

} 

f or(ctrl=l ;ctrl<9 ;ctrl++){ 

fscanf  (fdefmodel,"  '/,lf  "  ,4model[ctrl]  )  ; 

>  /*  for  Ctrl  */ 
fclose(fdefmodel) ; 

/*  model [1]  =  n,  nom  .99  */ 
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/*  model [2]  =  m,  nom  1.4  */ 

/♦  model [3]  =  mu,  nom  .002  */ 

/*  model [4]  =  sigma,  nom  .01  */ 

/*  model [5]  =  phi,  no  effect  for  delx=dely  */ 

/♦  model [6]  =  delx,  nom  .05236  radians  ♦/ 

/*  model [7]  =  dely,  nom  .05236  radians  */ 

/*  model [8]  =  k,  nom  5.0265  for  1200  Hz  Q  1500  m/s  */ 


/***♦**♦  MAIN  PROGRAM  LOOP  ♦♦♦♦♦♦♦/ 

/*  Print  net  arch  to  stdio  */ 
printf("Net  Architecture: \n") ; 
printf("No.  layers  *  V.ldXn"  ,arch[0]  )  ; 
printf  ("Layer  1  nodes  *  '/.ld\n"  ,arch[l]  )  ; 
printf  ("Layer  2  nodes  =  */,ld\n"  ,arch[2]  )  ; 
printf  ("Layer  3  nodes  =  ‘/,ld\n"  ,arch[3]  ) ; 
printf  ("No.  inputs  =  ‘/,ld\n\n",arch[4])  ; 

/*  Print  theta  range  to  stdio  ♦/ 
printf  ("Theta  paraimetersrXn") ; 
printf  ("min  =  '/,ld\n",  min); 
printf ("max  *  Xld\n",  max); 
printf  ("step  *  */,lf  \n\n"  ,  step); 

/*  Print  default  model  parameters  ♦/ 
printf  ("Default  model  pareLmeters:\n")  ; 
printf ("n  =  Xlf \n" ,model [1] ) ; 

printf  ("m  =  '/,lf\n",model[2])  ; 

printf ("mu  =  Xlf \n" ,model [3] ) ; 
printf  ("sigma  =  '/,lf\n"  ,model[4]  )  ; 
printf  ("del  =  '/,lf  \n\n"  ,model  [6]  )  ; 

/*  Print  the  number  of  training  sets  to  stdio  */ 
printf  ("N\imber  of  sets  =  ’/.IdXn"  ,num_sets)  ; 

f or(set_cnt=0 ; set_cnt<num_sets ; set_cnt++){ 

/*  This  is  loop  for  number  of  training  sets  */ 

/*  Get  new  training  set,  after  first  pass  only  */ 
if (set_cnt>0){ 

f or(ctrl=0 ;ctrl<5;ctrl++){ 
flag  »  fscanf(f train,"  */,lf " ,&train[ctrl] )  ; 

}  /*  for  Ctrl  */ 


if (flag  !=  1){ 

printf("EOF  reached  in  train.dat:  terminating\n") : 
return ; 

}  /*  if  flag  */ 

}  /♦  if  set_cnt  */ 

/♦  Print  training  setup  for  this  training  set  */ 
printf("Set  No.  ‘/.ld\n"  ,set_cnt)  ; 
printf  ("Learn  Gain  =  ‘/.If  \n"  .train [0]  )  ; 
printf  ("Learn  Switch  =  '/.If \n"  .train [1]  )  ; 
printf  ("Anneal  Gain  =  '/.If  \n"  .train[2]  )  ; 
printf  ("Init  Switch  =  '/.If  \n" .train[3]  )  ; 
printf  ("No.  reps  =  '/.If  \n"  .train[4]  )  ; 

/*  Assign  iteration  error  memory  */ 

iterr  =  (double  -J-)malloc(sizeof  (double)*(int)train[4]  )  ; 
if(iterr  ==  NULL){ 
printf ("iterr  allocation  errorXn") ; 
return; 

} 

/*  Initialize  error  */ 
perr  =  0; 

f or (t emp=0 ; t emp< (long) train [4] ; temp++) { 

/*  Loop  for  repetition  of  sajne  training  set  */ 

/*  Obtain  input  and  estimated  data  */ 

/♦  Compute  a  random  value  for  delx&dely,  scale  in  the  range 
of  .01745  to  .08727  for  the  desired  signal  (1  to  5  degrees)  */ 

do{ 

model  [6]  =  (( (double) rand () /RAND_MAX) ) ;  /*  range  is  0  to  1  */ 
model [6]  =  .08727  *  model  [6]; 

}while(model[6]  <  .01745); 

model [7]  =  model [6]; 
des[0]  =  model  [6]; 

/♦  Randomize  the  nm  parameters  to  try  and  make  del  estimation 
independent  of  them  */ 

do{ 
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model[2]  =  (((double)raiid()/RAND_MAX))  ;  /*  range  is  0  to  1  */ 
model [2]  =  .8  +  model [2]; 

}while(model [2]  <  1.2);  /*  range  is  1.2  to  1.8  */ 

model [l]  =  model [2]  ♦  .3833  +  .51;  /♦  range  is  .97  to  1.2  */ 

/♦  Randomize  the  mu  paraimeter:  range  .0002  tO  .02  */ 
do{ 

model[3]  =  ( ((double) rand () /RAND .MAX) ) ;  /♦  range  is  0  to  1  */ 
model [3]  =  . 02  *  model [3] ; 

}while(model[3]  <  .0002); 

/♦  Call  bism.c  to  generate  input  vector,  scale  for  0  to  1  */ 

/*  theta  is  from  71  to  88  degress  for  del  test,  89  degrees 
gives  too  large  a  value  of  ss  (2.5)  for  the  program  to  handle. 
Using  71  degrees  vice  60  to  reduce  sensitivty  to  changes  in  mu. 
Larger  errors  were  obtained  at  60  deg  even  when  noramlizing  the 
curve  to  the  start  angle  magnitude,  and  this  is  appraently  due 
to  the  curvature  caused  by  mu  ♦/ 

flag  *  bism(ss,  theta,  model); 
if (flag  1){ 

printf ("Error  in  bism  sub\n"); 
return ; 

} 

f or (ctrl=0 ; ctrKmodel [0]  ; ctrl++){ 

in[ctrl]=ss[ctrl]-ss[0] ;  /*  normalize  ampl  for  71  deg  */ 

} 

in[0]=l;  /*  Give  a  constant  to  play  with  */ 

/♦  Call  to  network  subroutine  */ 
nn(in,des ,est , arch, train, taps) ; 

/♦  printf  ("interset  '/.Id,  est  =  '/le\n"  ,temp,est  [0]  )  ;  */ 


/♦  Computer  errors  */ 

iterr[temp]  =  f abs (des [0] -est [0] )*100 .  /  aes [0] ; 
perr  +=  iterr[temp]; 


}  /*  for  temp  */ 

/*  Compute  err  std  deviation  */ 
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perr  =  perr  /  (double)train[4] ; 
errdev  =  0; 

for(ctrl=0;ctrl<train[4] ;ctrl++){ 

errdev  +=  (iterr[ctrl]  -  perr)  *  (iterr[ctrl]  -  perr); 

} 

errdev  =  sqrt(errdev  /  ((double)train[4]-l)) ; 

/♦  Print  results  of  network  for  end  of  current  set  */ 

printfC'avg  perr  =  ‘/,lf\n"  ,perr)  ; 
printf("perr  dev  =  '/.If \n"  .errdev)  ; 
printf("last  des  =  '/.If  \n"  ,des  [0]  )  ; 
printf("last  est  =  '/.lf\n"  .estCO])  ; 
printf("last  n  =  '/.If  \n"  .model  [l] ) ; 
printf("last  m  =  '/.If  \n"  .model  [2]  )  ; 
printf("last  mu  =  '/.If  \n"  .model  [3]  )  ; 

/*  Print  errors  to  ascii  file  ♦/ 
fiterr  =  f open("iterr. out" . "w") ; 
if(fiterr  ==  NULL){ 

printf ("Cannot  open  iterr  file  {w}\n"); 

}  /*  if  fiterr  =*  NULL  */ 
else{ 

f or(ctrl=0 ; ctrl<train[4] ; ctrl++)-( 
fprintf  (fiterr ."'/.If \n"  . iterr [ctrl] )  ; 

}  /*  for  Ctrl  */ 

}  /♦  else  ftaps_out  */ 
fclose(f iterr) ; 

free(iterr) ; 

}  /♦  for  set_cnt  */ 

/*^:i**^***  Finish  up  ★♦***♦***/ 

/*  Close  training  file  */ 

fclose(ftrain) ; 

/♦  Save  arch  and  taps  to  output  file  */ 
ftaps_out  =  f open("taps .out" ."w") ; 
if(ftaps_out  ==  NULL){ 
printf ("Cannot  open  taps  file  {w}\n"); 

}  /*  if  ftaps.out  ==  NULL  */ 
else{ 
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printf ("Dumping  arch  and  taps\n\n\n") ; 
for(ctrl=0; ctrl<5 ;ctrl++){ 
f printf  (f taps_out , "  ‘/.Id" ,  arch  [ctr l]  )  ; 

} 

f  or  (ctrl=0 :  ctrKno.layers ;  ctrl++)  { 

f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 
f or (ctr3=0 ; ctr3<no_inputs [ctrl] ; ctr3++) { 
f  printf  (ftaps.out,"  '/.le"  .taps  [ctrl]  [ctr2]  [ctr3] ) ; 
}  /♦  for  tap  */ 

}  /*  for  node  */ 

}  /♦  for  layer  */ 

}  /♦  else  ftaps.out  ♦/ 
f close(ftaps_out) ; 

/*  Free  memory  */ 
free(in) ; 
free(des) ; 
free(est) ; 
free(arch) ; 
free (train) ; 
free(theta)  ; 
free(ss) ; 

/*  free  tap  memory  */ 
f or (ctrl*0 ; ctrl<no_layers ; ctrl++) { 

f or (ctr2=0 ; ctr2<no_neurons [ctrl]  ; ctr2++) { 
free (taps [ctrl] [ctr2]  ) ; 

} 

} 

for(ctrl=0;ctrl<no_layers;ctrl++){ 
free(taps[ctrl]  ) ; 

} 

f ree(taps) ; 

return; 

} 
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Appendix  C 

BISSM  Subroutine 


/♦  Brian  Bourgeois,  Written  26AUG91,  Modified  26AUG91  */ 

/*  Program  bism.c  */ 

/♦  This  program  will  compute  the  various  parts,  and  the  total 

scattering  strength  based  upon  the  bissm2  model  and  its  various 
parameters.  The  output  is  a  vector  of  scattering  values  for 
a  specified  range  of  theta  (incident  ray  angle) .  */ 

#include  "bism.h" 

int  bisro(ss,  theta,  model) 
double  *ss, 

*theta, 

^fflodel ; 

{ 

/*  Model  Parameters 

model [0]  *  size  of  data  vectors  theta  and  ss 
model [l]  =  n.  Ratio  of  sound  speeds 

model [2]  =  m.  Ratio  of  densities 

model [3]  *  mu,  Lambert  coefficient 

model [4]  =  sigma.  Microscale  heights  roughness 
model [5]  =  phi,  azimuth  angle  in  radians 
model [6]  =  delx,  RMS  slope  deviation  along  track,  radiains 

model [7]  *  dely,  RMS  slope  deviation  across  track,  radians 

model [8]  *  k ;  Wavenumber 

*/ 

/*  Model  intermediate  variables  ♦/ 
double  tempi,  /*  interra  variable  */ 

temp2,  /*  interm  variable  */ 

temp3,  /♦  interm  variable  */ 

R,  /*  Rayleigh  coefficient  */ 
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R_a,  /*  part  of  Rayleight  computation  */ 

R_b,  /*  part  of  Rayleight  computation  ♦/ 

g.  /*  g  */ 

ae,  /*  angular  dependent  exponential  part  */ 

ref,  /*  reflective  component  */ 

lam;  /♦  Lambert  component  */ 

/♦  Misc  variables  */ 
long  cntl,  /*  counter  */ 
cnt2,  /*  counter  */ 
flag;  /*  error  flag  */ 

/♦  Note:  The  smallest  value  of  thet  must  be  large  enough  with 
respect  to  n  for  the  Rayleigh  reflection  coefficient  calculation. 
This  can  be  made  program  adjustable  in  the  future,  but  for  now  lets 
just  flag  it  */ 

if((modelCl]  <  1)  kk  (theta[0]  <  acos (model [1] ) )){ 
printf ("Error  in  theta  rangeXn") ; 
flag  =  1; 
return  flag; 

} 

/*  Scattering  Stength  Computations  */ 

f  or(cntl=0;  cntl<model  [O]  ;  cnt  !++)•( 

/*  Rayleigh  Coefficient  */ 

R_a  =  model [2]  ♦  sin(theta[cntl] ) ; 
tempi  =  model [l]  *  model [l]; 
temp2  =  cos (thetaCcntl] ) ; 
temp2  =  temp2  *  temp2; 

R_b  =  sqrt (tempi  -  temp2) ; 

R  *  (R_a  -  R.b)/(R.a  +  R_b) ; 

/*  compute  g  */ 

tempi  =  2  ♦  model[4]  *  modelCs]  *  sin(theta[cntl] ) ; 
g  =  tempi  *  tempi; 

/*  Compute  amgular  dependent  exp  part  */ 
tempi  =  cos (model [5] ) ; 

tempi  =  (tempi  *  templ)/(model[6]*model[6] ) ; 
temp2  *  sin(model [5] ) ; 

temp2  *  (temp2  *  temp2)/(model[7]*model[7] ) ; 
tempi  =  -.5  *  (tempi  +  temp2) ; 
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temp2  =  cos (thetaCcnt 1] )/sin(theta Cent  1] ) ; 
tempi  *  tempi  ♦  (temp2  *  temp2) ; 

tempi  *  exp(templ)/(8  *  PI  ♦  model [6]  *  model [7]); 

temp2  *  sin(theta[cntl] ) ; 

temp2  =  (temp2  *  temp2)  *  (temp2  *  temp2) ; 

ae  =  templ/temp2; 

/*  Compute  reflective  component  */ 

ref  =  (R  ♦  R)  ♦  exp(-l*g)  ♦  ae; 

/♦  Compute  lambert  component  */ 

tempi  =  s in  (theta  Cent  1]  )  ; 

leim  =  model  [3]  *  (tempi  *  tempi); 

/*  Compute  total  scattering  strength  ♦/ 
ssCcntl]  =  ref  +  lam; 

}  /*  for  cntl,  scatter  stength  loop  */ 

return; 

> 
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Appendix  D 

Parameter  Estimation  Programs 

D.l  fj.  Parameter  Estimation 


/♦  muest . c 

Brian  Bourgeois  Created;  25SEP91  Last  Mod:  25SEP91 

This  program  is  used  to  estimate  the  mu  parameter  of  the 
BISSM  algorithm  given  a  vector  that  is  scattering  strength 
vs.  angle  of  incidence.  The  program  loads  its  network 
architecture  and  taps  from  file  taps. in.  Note  that  the 
number  of  ss  points  provided  to  this  prograun,  and  their 
respective  amgles  is  fixed  by  the  network  architecture  and 
its  training.  An  area  is  provided  in  this  program  to  call 
for  data  input.  */ 

#define  RAND.MAX  2147483647  /*  2  31  -1  */ 

#define  PI  3.14159265358979 

#include  <stdio.h> 

#include  <malloc.h> 

#include  <math.h> 
tinclude  <ctype.h> 

#include  <stdlib.h> 

extern  int  nn (double*,  double*,  double  *,  long  *, 
double  *,  double  ***) ; 

extern  int  bism(double  *ss,  double  *theta,  double  *model) ; 

/*  exclude  this  for  cc  */ 

main (int  argc,  char  *argv[]) 

/*  main (argc, argv) 
int  argc; 


char  ♦argv [] ;  ♦/ 

{ 

/4i*4i4<***  Declare  variables  *♦♦♦*♦♦♦/ 


/*  network  configuration  */ 
long  arch [5];  /*  network  configuration  parameters  */ 

double  train [5] ;  /*  network  parameters  */ 
double  ***taps;  /*  network  weights  */ 

long  no_layers;  /*  Number  of  network  layers  with  neurons  */ 

long  *no_neurons;  /*  Number  of  neurons  in  each  layer  */ 

long  *no_inputs;  /*  Number  of  inputs  for  a  neuron  in  a  given  layer  */ 


/*  data  variables  */ 

double  des;  /*  network  desired  signal  */ 
double  est ;  /*  network  estimated  signal  ♦/ 

double  perr;  /*  desired  signal  squared  */ 
double  errdev;  /*  error  std  deviation  */ 
double  *iterr;  /*  Error  for  each  iteration  */ 


/♦  Misc  variables  */ 


long 


double 

FILE 

long 


double 


Ctrl,  /*  counter  ♦/ 

ctr2,  /*  counter  */ 

ctr3,  /*  counter  */ 

flag, 

temp ; 

ftemp; 

*ftaps_in;  /*  Input  taps  file,  with  arch  header  */ 

cntl,  /*  counter  */ 

min,  /*  minimum  thet  angle  */ 

max;  /*  maximum  thet  angle  */ 

step;  /*  thet  angle  step  size  */ 


/*  Model  variables  */ 


double  *theta,  /*  incident  angle  array  */ 
♦ss,  /♦  scattering  strength  data  */ 
model [9];  /*  model  parameters  ♦/ 

/*  Model  Paraimeters 

model [0]  =  size  of  data  vectors  theta  and  ss 
model [1]  =  n.  Ratio  of  sound  speeds 

model [2]  =  m.  Ratio  of  densities 

model [3]  =  mu,  Lambert  coefficient 
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model  [4]  =  sigma, 
model [5]  =  phi, 
model [6]  *  delx, 
model [7]  =  dely, 
model  [8]  =  k; 

*/ 


Microscale  heights  roughness 

azimuth  angle  in  radians 

PIMS  slope  deviation  along  track,  radians 

RMS  slope  deviation  across  track,  radians 

Wavenumber 


/*  Initialize  variables  */ 
flag  =  0; 

train C0]=0;  /*  Not  used  */ 
train Cl]=0;  /*  Not  used  */ 
train [2] =0;  /*  Not  used  */ 
train [3] =0;  /*  Not  used  */ 

train [4] =1000;  /*  No.  of  iterations  to  do  */ 

/*  Compute  theta  vector  length  */ 

/*  This  must  match  the  training  motif  used  for  the  network  in  use  */ 

min  =  15;  /*  minimum  theta  in  degrees  */ 

max  =60;  /♦  max  theta  */ 

step  =  2;  /♦  theta  angle  step  size  ♦/ 

/*  compute  vector  length  ♦/ 

model  Co]  =  (double) (max  -  min  +  l)/step  ;  /*  size  of  theta  vector  */ 

/*♦♦♦**  Obtain  network  configuration  data  ******/ 

/*  Read  in  architecture  from  taps  file,  taps. in  */ 
ftaps_in  =  f open("taps . in" , "r") ; 
if(ftaps_in  ==  NULL){ 

printf ("taps . in  does  not  exist  —  terminating\n") ; 
return; 

} 

for ( Ctrl =0; ctrl<5;ctrl++){ 

f  scanf  (ftaps_in,"  */,ld"  ,&arch[ctrl]  )  ; 

>  /*  for  Ctrl  ♦/ 

/*  archCO]  =  Number  of  layers  */ 

/*  archCl]  =  Number  of  neurons  in  first  layer  */ 

/*  arch [2]  =  Number  of  neurons  in  second  layer  */ 

/♦  arch [3]  =  Number  of  neurons  in  output  layer  */ 

/*  arch [4]  =  Number  of  inputs  to  first  layer  */ 

/♦  Load  up  architecture  variables  */ 

no.layers  =  arch[0]; 
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no_neurons  =  (long  *)malloc(sizeof (long)*(int)arch[0] ) ; 
if (no.neurons  ==  NULL){ 
printf ("no_neurons  allocation  error\n"); 
return; 

} 

no_neurons [0]  =  arch[l] ; 
no_neurons[l]  =  arch [2] ; 
no.neurons [2]  =  arch [3] ; 

no_inputs  =  (long  *)malloc(sizeof (long)* (int) arch [O] ) ; 
if(no_inputs  ==  NULL){ 
printf ("no_inputs  allocation  errorXn"); 
return; 

} 

no_inputs[0]  =  arch [4] ;  /*  Inputs  to  first  layer  */ 

/♦  Note  this  should  be  same  as  length  of  ss  vector  */ 
no_inputs[l]  =  no_neurons [0] ;  /*  Inputs  to  second  layer  */ 

/*  same  as  no_neurons  in  first  layer  for  full  connectivity  */ 
no_inputs[2]  =  no_neurons  [1] ;  /*  Inputs  to  output  layer  */ 

/*  saime  as  no^neurons  in  second  layer  for  full  connectivity  */ 

/******  Allocate  storage  ******/ 

theta  =  (double  *)malloc(sizeof (double)*(int)modelCO] ) ; 
if (theta  =*  NULL){ 

printf ("theta  allocation  errorXn"); 
return; 

} 

ss  =  (double  *)malloc(sizeof (double)*(int)model [O] ) ; 
if(ss  ==  NULL){ 

printf ("ss  allocation  errorXn") ; 
return; 

} 

iterr  =  (double  *)malloc(sizeof (double)*(int)train[4] ) ; 
if(iterr  ==  NULL){ 
printf ("iterr  allocation  errorXn") ; 
return; 

} 

/*  nn  weights  allocation  is  a  bit  more  involved  */ 
taps  =  (double  ♦♦*)malloc(sizeof (double  ***)*no_layers) ; 
if (taps  ==  NULL){ 
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printf("taps  allocation  error,  level  l\n") ; 
return; 

}  /♦  if  taps  ==  NULL  */ 
f or (ctrl=0 ; ctrl<no_layers ; ctrl++) { 

taps[ctrl]  =  (double  **)raalloc(si2eof (double  ♦*) *no_neurons [ctrl] ) ; 
if (taps [ctrl]  ==  NULL){ 

printf("taps  allocation  error,  level  2\n"); 
return; 

}  /♦  if  taps[]  ==  NULL  */ 

}  /*  for  Ctrl  */ 

for(ctrl=0;  ctrl<no_layers  ;ctrl++)-C 
f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 

taps [ctrl] [ctr2]  =  (double  ♦)malloc(sizeof (double) *no_inputs [ctrl] ) 
if (taps [ctrl] [ctr2]  ==  NULL){ 
printfC'taps  allocation  error,  level  3\n"); 
return; 

}•  /♦  if  taps[][]  ==  NULL  */ 

}  /♦  for  ctr2  */ 

}  /*  for  Ctrl  */ 

/>lcit!**it*****t  ******************  ******************/ 

/*  Working  Area  of  MAIN  */ 

/******  Load  taps  array  *****=i</ 

/*  Load  taps  from  taps. in  and  close  file  */ 
f  or  (ctr  1=0;  ctrl<no_  layers ;  Ctrl ++)■( 

f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 
f or (ctr3=0;ctr3<no_ inputs [ctrl] ;ctr3++){ 
fscanf  (ftaps.in,"  */,le"  ,&taps  [ctrl]  [ctr2]  [ctr3]  ) ; 

}  /*  for  tap  */ 

)■  /*  for  node  */ 

}  /*  for  layer  */ 
fclose(ftaps_in) ; 

/*  Assign  incident  angles  */ 
for(cntl=0; cntl<model [0] ; cntl++){ 

theta[cntl]  =  (double) (min  +  cntl*step)  *  PI  /  180.; 

/*  printf  ("theta  =  '/,f  \n"  ,  theta [cntl]  )  ;  */ 

> 

/*  Load  nominal  model  parameters  */ 


model  [l]  =  .99; 


model [2]  =  1.4; 
model [3]  =  .002; 
model [4]  =  .01; 
model [5]  =  0 ; 
model [6]  =  .05236; 
model [7]  =  .05236; 
model [8]  =  5.0265; 

/**^f****  main  program  loop  ♦♦♦♦♦♦♦/ 

/*  Loop  for  multiple  estimations  ♦/ 
printf ("WORKINGXn") ; 

f or(temp=0;temp< (long) train [4] ;temp++) { 

/*  Obtain  input  and  estimated  data  ♦/ 

/*  Compute  a  random  value  for  delx&dely,  scale  in  the  range 
of  .01745  to  .08727  for  the  desired  s..^nal  (l  to  5  degrees)  */ 

do{ 

model[6]  =  ( ((double) rand () /RAND_MAX) ) ;  /*  range  is  0  to  1  */ 
model [6]  *  .08727  ♦  model  [6] ; 

)-while(model  [6]  <  .01745); 

model  [7]  =  model  [6]; 

/♦  Compute  a  rcindora  value  for  nra.  Note  that  n  amd  m  are 
computed  as  being  linearly  related  */ 

do{ 

model [2]  =  ( ((double) rand () /RAND.MAX) ) ;  /*  range  is  0  to  1  */ 
model [2]  =  .8  +  model  [2]; 

>while(model [2]  <  1.2);  /*  range  is  1.2  to  1.8  */ 

model [1]  =  model [2]  ♦  .3833  +  .51;  /*  range  is  .97  to  1.2  */ 

/*  Randomize  the  mu  parameter:  range  .0002  tO  .02  ♦/ 
do{ 

model[3]  =  ( ((double) rand () /RAND_MAX) ) ;  /♦  range  is  0  to  1  */ 
model [3]  =  .02  *  model  [3]; 

}while(model [3]  <  .0002); 

des  =  model [3] ; 
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/*  Call  bism.c  to  generate  input  vector,  scale  for  range  0  to  1  */ 

flag  =  bism(ss,  theta,  model); 
if (flag  ==  1){ 

printf ("Error  in  bism  sub\n") ; 
return ; 

} 

/♦  Adjust  ss  data  for  proper  network  operation,  as  per  training  */ 
/*  For  10  to  60  deg,  ss  rainges  from  3*10“-6  to  3*10“-2  */ 

/*  Scale  up  from  .03  to  .3  */ 
f or(ctrl=0; ctrl<model [0] ; ctrl++){ 
ss[ctrl]=ss[ctrl]*10.  -  .15; 

} 

/♦  Call  to  network  subroutine  */ 


nn(ss ,&des,&est , arch, train, taps) ; 

/*  Compute  errors  */ 

iterr[temp]  =  fabs(des-est)*100./des; 

}  /*  for  temp  ♦/ 

/*  Compute  err  and  its  std  deviation  */ 
perr  =  0; 

f or (ctrl=0 ; ctrl<train[4] ; ctrl++){ 
perr  +=  iterr[ctrl]; 

} 

perr  =  perr  /  (double)train[4] ; 
errdev  =  0 ; 

f or (ctrl=0 ; ctrl<train[4] ; ctrl++){ 

errdev  +=  (iterr[ctrl]  -  perr)  *  (iterr[ctrl]  -  perr); 

} 

errdev  =  sqrt(errdev  /  ((double)train[4]-l)) ; 


/♦  Print  results  of  network  test  ♦/ 


printf ("avg  perr 
printf ("perr  dev 
printf ("last  des 
printf ("last  est 


'/.If \n", perr)  ; 
'/,lf\n" ,  errdev) ; 
'/.lf\n"  ,des)  ; 
'/,lf\n"  ,est)  ; 


/********  Finish  up  *********/ 


53 


/♦  Free  memory  ♦/ 
free(iterr) ; 
free (theta) ; 
free(ss) ; 
free(no_neurons) ; 
free(no_inputs) ; 

/*  free  tap  memory  */ 
f or (ctrl~0 ; ctrl<no_layers ; ctrl++) { 

f or (ctr2=0 ; ctr2<no .neurons [ctrl] ; ctr2++) { 
free  (taps  [ctrl]  [ctr2]  )  ; 

} 

} 

f or (ctrl=0 ; ctr l<no_layers ; ctrl++) { 
free (taps [ctrl] ) ; 

} 

free (taps) ; 

return ; 

} 


D.2  nm  Parameter  Estimation 

/*  nmest.c 

Brian  Bourgeois  Created:  25SEP91  Last  Mod:  25SEP91 

This  progreim  is  used  to  estimate  the  nm  parameter  of  the  BISSM 
algorithm  given  a  vector  that  is  scattering  strength  vs.  angle 
of  incidence.  The  program  loads  its  network  architecture  and 
taps  from  file  taps. in.  Note  that  the  number  of  ss  points 
provided  to  this  program,  and  their  respective  angles  is  fixed 
by  the  network  architecture  and  its  training. 

An  area  is  provided  in  this  program  to  call  for  data  input. 

*/ 

#define  RAND.MAX  2147483647  /*  2  31  -1  */ 

#define  PI  3.14159265358979 

#include  <stdio.h> 

♦include  <malloc.h> 

♦include  <math.h> 

♦include  <ctype.h> 
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#include  <stdlib.h> 

extern  int  nn (double*,  double*,  double  *,  long  *, 
double  *,  double  ***) ; 

extern  int  bism(double  *ss,  double  *theta,  double  *model) ; 

/*  exclude  this  for  cc  */ 

mainCint  argc,  char  *argvC]) 

/*  main ( argc, argv) 

int  argc; 

char  *argv[];  */ 

{ 

/*******  Declare  variables  *♦♦♦**♦*/ 

/*  network  configuration  */ 
long  arch [5];  /*  network  configuration  parameters  */ 

double  train [5];  /*  network  parameters  */ 
double  ***taps;  /*  network  weights  */ 

long  no.layers;  /*  Number  of  network  layers  with  neurons  */ 

long  *no_neurons;  /*  Number  of  neurons  in  each  layer  */ 

long  *no_inputs;  /*  Number  of  inputs  for  a  neuron  in  a  given  layer  */ 

/*  data  variables  */ 

double  des;  /*  network  desired  signal  */ 
double  est;  /*  network  estimated  signal  */ 
double  perr;  /*  desired  signal  squared  */ 
double  errdev;  /*  error  std  deviation  */ 
double  *iterr;  /*  Error  for  each  iteration  */ 

/*  Misc  variables  */ 

long  Ctrl,  /*  counter  */ 

ctr2,  /*  counter  */ 
ctr3,  /*  counter  */ 
flag, 
temp; 

double  ftemp; 

FILE  *ftaps_in;  /*  Input  taps  file,  with  arch  header  */ 
long  cntl,  /*  counter  */ 

min,  /*  minimum  thet  angle  */ 

max;  /*  maximum  thet  angle  */ 

double  step;  /*  thet  angle  step  size  */ 

/♦  Model  variables  */ 
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double  ♦theta,  /♦  incident  angle  array  */ 
*ss,  /*  scattering  strength  data  ♦/ 
model [9];  /*  model  parameters  ♦/ 


/♦  Model 
model [0] 
model [1] 
model [2] 
model [3] 
model [4] 
model [5] 
model [6] 
model [7] 
model [8] 
*/ 


Parauneters 

=  size  of  data  vectors  theta  and  ss 
=  n,  Ratio  of  sound  speeds 

=  m,  Ratio  of  densities 

=  mu,  Lambert  coefficient 
=  sigma.  Microscale  heights  roughness 
=  phi,  azimuth  angle  in  radians 
=  delx,  RMS  slope  deviation  along  track,  radians 
=  dely,  RMS  slope  deviation  across  track,  radians 

=  k;  Wavenumber 


/♦  Initialize  variables  */ 
flag  =  0; 

train[0]=0;  /♦  Not  used  ♦/ 

train[l]=0;  /*  Not  used  ♦/ 

train[2]=0;  /♦  Not  used  ♦/ 

train [33=0;  /♦  Not  used  */ 

train[4]*1000;  /♦  No.  of  iterations  to  do  */ 

/*  Compute  theta  vector  length  */ 

/♦  This  must  match  the  training  motif  used  for  the  network  in  use  */ 

min  =71;  /*  minimum  theta  in  degrees  */ 

max  =88;  /♦  max  theta  */ 

step  =  1;  /*  theta  angle  step  size  */ 

/♦  compute  vector  length  ♦/ 

model[0]  =  (double) (max  -  min  +  l)/step  ;  /*  size  of  theta  vector  ♦/ 

/******  Obtain  network  configuration  data  *♦**♦*/ 

/♦  Read  in  architecture  from  taps  file,  taps. in  */ 
ftaps_in  =  f open("taps . in" , "r") ; 
if(ftaps_in  ==  NULL){ 

printf ("taps . in  does  not  exist  —  terminatingXn") ; 
return; 

> 

f or(ctrl=0; ctrl<5 ;ctrl++){ 

f  sc2uif  (ftaps_in, "  '/.Id"  ,&arch[ctrl]  )  ; 

}  /*  for  Ctrl  ♦/ 
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/♦  arch[0]  =  Number  of  layers  */ 

/♦  archCl]  =  Number  of  neurons  in  first  layer  */ 

/♦  arch [2]  =  Number  of  neurons  in  second  layer  */ 

/♦  arch [3]  =  Number  of  neurons  in  output  layer  */ 

/♦  arch [4]  =  Number  of  inputs  to  first  layer  */ 

/♦  Load  up  architecture  variables  */ 

no.layers  =  arch[0] ; 

no_neurons  =  (long  *)malloc(sizeof (long)*(int)arch[0] ) ; 
if (no.neurons  ==  NULL){ 
printfC "no .neurons  allocation  errorXn"); 
return; 

} 

no.neurons [0]  =  archCl] ; 
no_neurons[l]  =  arch [2] ; 
no.neurons [2]  =  arch [3] ; 

no. inputs  =  (long  ♦)malloc(sizeof (long)*(int)arch[0] ) ; 
if (no. inputs  ==  NULL){ 
printf ("no. inputs  allocation  errorXn"); 
return; 

} 

no.inputsCO]  *  arch [4] ;  /*  Inputs  to  first  layer  */ 

/*  Note  this  should  be  seune  as  length  of  ss  vector  */ 
no.inputsCl]  *  no.neurons [O] ;  /♦  Inputs  to  second  layer  ♦/ 

/♦  same  as  no.neurons  in  first  layer  for  full  connectivity  ♦/ 
no. inputs [2]  *  no.neurons [l] ;  /*  Inputs  to  output  layer  */ 

/*  same  as  no.neurons  in  second  layer  for  full  connectivity  * 

Allocate  storage  ♦****♦/ 

theta  =  (double  *)malloc(sizeof (double) *(int)model [0] ) ; 

if (theta  **  NULL){ 

printf ("theta  allocation  errorXn") ; 
return; 

} 

ss  =  (double  *)malloc(sizeof (double)*(int)model[0] ) ; 
if(ss  =*  NULL){ 

printf ("ss  allocation  errorXn"); 
return; 

} 


iterr  =  (double  *)malloc(sizeof (double)*(int)trainC4] ) ; 
if  (it err  ==  NULL)-C 
printf ("iterr  allocation  error\n") ; 
return ; 

} 

/*  nn  weights  allocation  is  a  bit  more  involved  ♦/ 
taps  =  (double  ***)malloc(sizeof (double  ♦**) *no_layers) ; 
if (taps  ==  NULL){ 

printf ("taps  allocation  error,  level  l\n") ; 
return; 

}  /*  if  taps  ==  NULL  */ 
f or (ctrl=0 ; ctr l<no_layers ; ctrl++) { 

tapsCctrl]  =  (double  **)malloc(sizeof (double  **) ♦no_neurons [ctrl] ) ; 
if (tapsCctrl]  ==  NULL){ 

printf ("taps  allocation  error,  level  2\n") ; 
return; 

}  /♦  if  taps[]  ==  NULL  */ 

}  /*  for  Ctrl  */ 

for(ctrl=0;ctrl<no_layers ;ctrl++){ 
f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 

taps [ctrl] [ctr2]  =  (double  *)malloc(sizeof (double) *no_inputs [ctrl]) ; 
if (taps [ctrl] [ctr2]  =*  NULL){ 
printf ("taps  allocation  error,  level  3\n"); 
return; 

}  /♦  if  taps[][]  ==  NULL  */ 

}  /*  for  ctr2  */ 

}  /*  for  Ctrl  */ 

!  *)([♦♦♦♦♦♦***♦♦****♦****  ******♦*♦*♦♦♦♦*  ***♦*♦*♦**/ 

/♦  Working  Area  of  MAIN  */ 

1 Load  taps  array  ♦***♦♦/ 

/♦  Load  taps  from  taps,  in  eoid  close  file  */ 
for (ctr 1=0; ctrl<no_layers;ctrl++){ 

f or (ctr2=0 ; ctr2<no .neurons [ctrl] ; ctr2++) { 
f or(ctr3=0 ; ctr3<no_inputs [ctrl] ;ctr3++) { 
f  scanf  (ftaps.in, "  ‘/,le"  ,&taps  [ctrl]  [ctr2]  [ctr3]  )  ; 

}  /*  for  tap 
>  /*  for  node  */ 

}  /*  for  layer  *! 
fclose(ftaps_in) ; 
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/*  Assign  incident  angles  ♦/ 
f or  (cnt  1=0 ;  cnt  l<niodel  [0]  ;  cnt  1++) { 

theta[cntl]  =  (double) (min  +  cntl^step)  *  PI  /  180.; 
/♦  printf  ("theta  =  */,f\n"  .theta [cnt  1]  )  ;  */ 

} 


/♦  Load  nominal  model  parameters  ♦/ 


model [1] 
model [2] 
model [3] 
model [4] 
model [5] 
model [6] 
model [7] 
model [8] 


.99; 

1.4; 

.  002 ; 
.01; 

0; 

.05236 

.05236 

5.0265 


/♦♦♦♦***  MAIN  PROGRAM  LOOP  ****♦♦*/ 


/♦  Loop  for  multiple  estimations  */ 
printf ("W0RKING\n") ; 

f or(temp=0;temp<(long)train[4] ;temp++){ 

/*  Obtain  input  and  estimated  data  ♦/ 

/*  Randomize  delx&dely,  scale  in  the  range  of  .01745  to  .08727 
(1  to  5  degrees)  */ 

do{ 

model [6]  =  ( ((double) rand () /RAND_MAX) ) ;  /*  remge  is  0  to  1  */ 
model [6]  =  .08727  ♦  model [6]; 

}while(model[6]  <  .01745); 

model [7]  =  model  [6]; 


/♦  Compute  a  raoidom  value  for  nm.  Note  that  n  and  m  are 
computed  as  being  linearly  related  */ 

do{ 

model[2]  =  (( (double) reuad () /RAND.MAX) ) ;  /*  range  is  0  to  1  */ 
model [2]  *  .8  +  model [2]; 

}while(model [2]  <  1.2);  /*  range  is  1.2  to  1.8  */ 
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model [l]  =  model [2]  *  .3833  +  .51;  /*  range  is  .97  to  1.2  */ 
des  *  model [2] /4.; 

/*  Randomize  the  mu  parameter:  raoige  .0002  tO  .02  */ 
do{ 

model[3]  =  ( ((double) rand () /RAND.MAX) ) ;  /♦  range  is  0  to  1  */ 
model [3]  =  .02  ♦  model  [3]; 

}while(modelC3]  <  .0002); 

/*  Call  bism.c  to  generate  input  vector,  scale  for  range  0  to  1  */ 

flag  =  bism(ss,  theta,  model); 
if  (flag  ==  1)-C 

printf ("Error  in  bism  sub\n"); 
return ; 

} 

/*  normalize  ampl  for  starting  incident  cingle  to  reduce 
mu  effect  on  estimate  */ 
ftemp  =  ss  [O]  ; 

f or (ctrl=0 ; ctrl<raodel [0] ; ctrl++){ 
ss [ctrl] ®ss [Ctrl] -ftemp ; 

} 

ss[0]=l;  /♦  Constant  input  node  for  network  */ 

/*  Call  to  network  subroutine  ♦/ 

nn(ss ,&des ,&est , arch, train, taps) ; 

/*  Compute  errors  */ 

iterrCtemp]  ®  f abs(des-est)*100 ./des; 

}  /*  for  temp  */ 

/♦  Compute  err  and  its  std  deviation  */ 
perr  *  0; 

for(ctrl=0;ctrl<train[4]  ;ctr !++)•( 
perr  +=  iterr[ctrl]; 

} 

perr  *  perr  /  (double) train [4] ; 
errdev  =  0; 

f or (ctrl*0 ; ctrl<train[4] ; ctrl++){ 

errdev  +*  (iterrCctrl]  -  perr)  *  (iterrCctrl]  -  perr); 

} 
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errdev  =  sqrt(errdev  /  ((double)train[4] -1)) ; 

/♦  Print  results  of  network  test  */ 

printf("avg  perr  =  */,lf  \n"  ,perr)  ; 
printf("perr  dev  =  */,lf\n"  , errdev)  ; 
printf("last  des  =  '/.lf\n"  ,4.  *des)  ; 
printf("last  est  =  */,lf\n"  ,4.*est)  ; 

Finish  up  *♦♦****♦♦/ 

/♦  Free  memory  */ 
free(iterr) ; 
free (theta) ; 
free(ss) ; 
f ree(no_neurons) ; 
f ree(no_inputs) ; 

/*  free  tap  memory  */ 
for(ctrl=0; ctrl<no_layers ;ctrl++){ 

f or(ctr2=0 ;ctr2<no_neurons [ctrl] ;ctr2++){ 
free  (taps  [ctrl]  [ctr2]  )  ; 

} 

} 

f or (ctrl=0 ; ctr l<no_layers ; ctrl++) { 
free (taps [ctrl] ) ; 

> 

f ree(taps) ; 

return; 

} 


D.3  S  Parameter  Estimation 

/*  delest.c 

Brian  Bourgeois  Created:  23SEP91  Last  Mod:  24SEP91 

This  program  is  used  to  estimate  the  del  parameter  of  the 
BISSM  algorithm  given  a  vector  that  is  scattering  strength 
vs.  angle  of  incidence.  Due  to  the  dependence  of  the  delx 
eoid  dely  parameters,  only  a  single  parameter,  del,  is  used 
wherein  delx  *  dely  *  del.  The  program  loads  its  network 
architecture  auid  taps  from  file  taps. in.  Note  that  the 
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number  of  ss  points  provided  to  this  progreun,  and  their 
respective  angles  is  fixed  by  the  network  architecture 
and  its  training.  An  area  is  provided  in  this  prograon  to 
call  for  data  input .  ♦/ 

#define  RAND.MAX  2147483647  /♦  2  31  -1  */ 

#define  PI  3.14159265358979 

♦include  <stdio.h> 

♦include  <malloc.h> 

♦include  <math.h> 

♦include  <ctype.h> 

♦include  <stdlib.h> 

extern  int  nn (double*,  double*,  double  *,  long  *, 
double  *,  double  ***) ; 

extern  int  bism(double  *ss,  double  *theta,  double  *model) ; 

/*  exclude  this  for  cc  */ 

mainCint  argc,  char  *argv[]) 

/*  main (argc, argv) 

int  argc; 

char  *argv[];  */ 

{ 

Declare  variables  *♦****♦*/ 

/*  network  configuration  */ 
long  arch [5] ;  /*  network  configuration  parameters  */ 

double  train [5] ;  /*  network  parameters  */ 
double  ***taps;  /*  network  weights  */ 

long  no.layers;  /*  Number  of  network  layers  with  neurons  */ 

long  *no_neurons;  /*  Number  of  neurons  in  each  layer  */ 

long  *no_inputs;  /*  Number  of  inputs  for  a  neuron  in  a  given  layer  */ 

/*  data  variables  */ 

double  des;  /*  network  desired  signal  */ 
double  est ;  /*  network  estimated  signal  */ 

double  perr;  /*  desired  signal  squared  */ 
double  errdev;  /♦  error  std  deviation  */ 
double  *iterr;  /*  Error  for  each  iteration  */ 

/*  Misc  variables  */ 

long  Ctrl,  /*  counter  */ 
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ctr2,  /*  counter  */ 
ctr3,  /*  counter  */ 
flag, 
temp ; 

double  ftemp; 

FILE  *ftaps_inj  /♦  Input  taps  file,  with  arch  header  */ 
long  cntl,  /*  counter  */ 

min,  /*  minimum  thet  angle  */ 

max;  /*  maximum  thet  angle  */ 

double  step;  /*  thet  angle  step  size  */ 

/*  Model  variables  */ 

double  ♦theta,  /*  incident  angle  array  */ 

*ss,  /*  scattering  strength  data  */ 
model [9];  /*  model  parameters  */ 

/*  Model  Parameters 

model  [0]  =  size  of  data  vectors  theta  and  ss 
model  Cl]  =  n.  Ratio  of  sound  speeds 

model [2]  =  m.  Ratio  of  densities 

model [3]  =  mu,  Lambert  coefficient 

model [4]  ®  sigma.  Microscale  heights  roughness 

model [5]  *  phi,  azimuth  angle  in  radians 

model [6]  =  delx,  RMS  slope  deviation  along  track,  rad*ans 

model [7]  =  dely,  RMS  slope  deviation  across  track,  radians 

model [8]  =  k;  Wavenumber 

*/ 

/*  Initialize  variables  */ 
flag  =  0; 

train [0]=0;  /*  Not  used  */ 

train[l]=0;  /*  Not  used  */ 

train[2]=0;  /*  Not  used  */ 

train [3] =0;  /*  Not  used  */ 

train[4]=1000;  /*  No.  of  iterations  to  do  ♦/ 

/♦  Compute  theta  vector  length  ♦/ 

/♦  This  must  match  the  training  motif  used  for  the  network  in  use  ♦/ 

min  =  71;  /*  minimum  theta  in  degrees  */ 

max  =88;  /*  max  theta  */ 

step  =  .5;  /*  theta  angle  step  size  */ 

/*  compute  vector  length  */ 

modelCO]  =  (double) (max  -  min  +  l)/step  ;  /*  size  of  theta  vector  */ 
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/**♦**♦  Obtain  network  configuration  data  ***♦♦*/ 

/*  Read  in  architecture  from  taps  file,  taps. in  */ 
ftaps_in  =  f  openC'taps  .  in"  .  "r")  ; 
if(ftaps_in  ==  NULL){ 

printf ("taps . in  does  not  exist  —  terminatingXn") ; 
return; 

} 

f or(ctrl=0 ; ctrl<5 ;ctrl++) { 

f  scanf  (ftaps.in, "  ‘/.Id"  ,&arch[ctrl]  ) ; 

}  /*  for  Ctrl  */ 

/*  arch[0]  =  Number  of  layers  */ 

/*  arch[l]  =  Number  of  neurons  in  first  layer  */ 

/*  arch [2]  =  Number  of  neurons  in  second  layer  */ 

/*  arch [3]  =  Number  of  neurons  in  output  layer  */ 

/*  archC4]  =  Number  of  inputs  to  first  layer  */ 

/*  Load  up  architecture  variables  */ 

no_layers  =  arch[0] ; 

no.neurons  =  (long  ♦)malloc(sizeof (long)* (int)arch [0] ) ; 
if (no.neurons  ==  NULL){ 

printf ("no.neurons  allocation  errorXn"); 
return; 

} 

no.neuronsCO]  =  arch[l] ; 
no_neurons [1]  =  arch [2] ; 
no_neurons [2]  =  arch [3] ; 

no.inputs  =  (long  *)malloc(sizeof (long)*(int)arch [0] ) ; 
if(no_i^puts  ==  NULL)-C 
printf ("no.inputs  allocation  errorXn") ; 
return; 

} 

no_inputsC0]  =  arch[4];  /*  Inputs  to  first  layer  */ 

/*  Note  this  should  be  same  as  length  of  ss  vector  */ 
no_inputs[l]  =  no.neurons [0] ;  /*  Inputs  to  second  layer  */ 

/*  same  as  no.neurons  in  first  layer  for  full  connectivity  */ 
no_inputs[2]  =  no_neurons [l] ;  /*  Inputs  to  output  layer  */ 

/*  same  as  no_neurons  in  second  layer  for  full  connectivity  */ 

/♦**♦**  Allocate  storage  ******/ 
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theta  =  (double  *)malloc(si2eof (double)*(int)model [0] ) ; 
if (theta  ==  NULL){ 

printf ("theta  allocation  errorXn") ; 
return; 

} 


ss  =  (double  *)malloc(sizeof (double)*(int)model [0] ) ; 
if(ss  »=  NULL){ 

printf ("ss  allocation  errorXn") ; 
return; 

} 

iterr  =  (double  *)malloc(si2eof (double) *(int)train [4] ) ; 
if(iterr  ==  NULL){ 
printf ("iterr  allocation  errorXn") ; 
return; 

} 

/♦  nn  weights  allocation  is  a  bit  more  involved  *! 
taps  =  (double  ***)malloc(si2eof (double  ***) *no_layers) ; 
if  (taps  =*  NULLX 

printf ("taps  allocation  error,  level  l\n"); 
return; 

}  /*  if  taps  ==  NULL  */ 
f or (ctrl=0 ; ctrl<no_layers ; ctrl++){ 

tapsLctrl]  =  (double  ♦♦)malloc(sizeof (double  **) *no_neurons [ctrl] ) ; 
if  (taps  [ctrl]  ==  NULD-C 

printf ("taps  allocation  error,  level  2\n") ; 
return; 

}  /*  if  taps[]  ==  NULL  */ 

}  /*  for  Ctrl  */ 

f or (ctrl=0 ; ctrl<no_ layers ;ctrl++) { 
f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 

taps [ctrl] [ctr2]  =  (double  *)malloc(si2eof (double) *no_inputs [ctrl] ) 
if (taps [ctrl] [ctr2]  ==  NULL){ 
printf ("taps  allocation  error,  level  3\n") ; 
return; 

}  /*  if  taps[][]  ==  NULL  */ 

>  /♦  for  ctr2  */ 

}  /♦  for  Ctrl  */ 

/♦*♦*****♦♦***♦******♦*****♦***♦****♦********+*♦/ 
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/♦  Working  Area  of  MAIN  */ 


Load  taps  array  ******/ 

/*  Load  taps  from  taps,  in  aind  close  file  */ 
f or(ctrl=0  ctrl<no_layers ;ctrl++){ 

f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) { 
f or (ctr3=0; ctr3<no_ inputs [ctrl] ;ctr3++){ 
fscajif  (ftaps_in,"  '/,le"  ,&taps  [ctrl]  [ctr2]  [ctr3]  )  ; 

}  /*  for  tap  */ 
y  /*  for  node  */ 
y  /*  for  layer  */ 
f close (ftaps_in) ; 

/♦  Assign  incident  angles  */ 
f or(cntl=0 ; cntl<model [0] ; cntl++){ 

thetaCcntl]  =  (double) (min  +  cntl*step)  *  PI  /  180.; 
/*  printf  ("theta  =  */,f\n"  , thetaCcntl]  )  ;  */ 

y 


/*  Load  nominal  model  parameters  ♦/ 


model [1] 
model [2] 
model [3] 
model [4] 
model [5] 
model [6] 
model [7] 
model [8] 


.99; 

1.4; 

.002; 

.01; 

0; 

.05236 

.05236 

5.0265 


Z^.******  main  program  loop  *♦****♦/ 


/*  Loop  for  multiple  estimations  */ 
printf ("WORKINGNn") ; 

for(temp=0 ;  temp<(long)train[4]  ;temp++)-C 
/*  Obtain  input  and  estimated  data  */ 

/♦  Compute  a  raindom  value  for  delxidely,  scale  in  the  range  of 
.01745  to  .0872  for  the  desired  signal  (1  to  5  degrees)  */ 

do{ 

model[6]  =  (( (double) rand() /RAND.MAX) ) ;  /*  range  is  0  to  1  */ 
model [6]  =  .08727  *  model  [6]; 
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}while(model [6]  <  .01745); 


model [7]  =  model [6]; 
des  =  model  [6] ; 

/*  Compute  a  random  value  for  nm.  Note  that  n  and  m  are 
computed  as  being  linearly  related  */ 

do{ 

model[2]  =  ( ((double) rand () /RAND.MAX) ) ;  /♦  ramge  is  0  to  1  */ 
model [2]  =  .8  +  model [2]; 

}while(model[2]  <  1.2);  /*  range  is  1.2  to  1.8  */ 

model  [l]  =  model  [2]  *  .3833  +  .51;  /♦  range  is  .97  to  1.2  ♦/ 

/*  Randomize  the  mu  narameter:  range  .0002  tO  .02  */ 
do{ 

model [3]  =  (( (double) rand ()/RAND_MAX)) ;  /*  range  is  0  to  1  */ 
model [3]  =  .02  *  model  [3]; 

}while(model  [3]  <  .0002); 

/♦  Call  bism.c  to  generate  input  vector,  scale  for  range  0  to  1  */ 

flag  =  bism(ss,  theta,  model); 
if (flag  ==  1){ 

printf ("Error  in  bism  sub\n"); 
return; 

} 

/*  normalize  ampl  for  starting  incident  Euigle  to  reduce 
mu  effect  on  estimate  ♦/ 
f temp  =  ss  [0]  ; 

f or(ctrl=0 ; ctrl<model [O] ; ctrl++){ 
ss [ctrl] =ss [ctrl] -f temp ; 

} 

ssC0]=l;  /*  Constant  input  node  for  network  */ 

/*  Call  to  network  subroutine  */ 

nn(ss,ftdes,ftest , arch, train, taps) ; 

/♦  Compute  errors  */ 

iterrCtemp]  =  f abs(des-est)*100 ./des ; 
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}  /♦  for  temp  */ 


/*  Compute  err  and  its  std  deviation  */ 
perr  =  0; 

f  or  (ctrl=0 ;  ctrl<train[4]  ;  ctrl++)-C 
perr  +=  iterr[ctrl] ; 

} 

perr  =  perr  /  (double) train [4] ; 
errdev  ®  0^ 

f or(ctrl=0;ctrl<train[4] ; ctrl++){ 

errdev  +=  (iterr[ctrl]  -  perr)  ♦  (iterr[ctrl]  -  perr); 

> 

errdev  =  sqrt(errdev  /  ( (double)train[4] -1) ) ; 

/♦  Print  results  of  network  test  */ 

printf("avg  perr  =  */,lf  \n"  ,perr)  ; 
printf("perr  dev  =  y,lf\n"  .errdev)  ; 
printf("last  des  =  */,lf  \n"  ,des)  ; 
printf("last  est  *  ‘/.If  \n"  ,est)  ; 

/♦♦♦♦♦♦♦♦  Finish  up  ♦*♦*♦♦***/ 

/*  Free  memory  */ 
f ree(iterr) ; 
free (theta) ; 
f ree(ss) ; 
f ree(no_neurons) ; 
free(no_inputs) ; 

/*  free  tap  memory  */ 
f or(ctrl=0;ctrl<no_layers ;ctrl++){ 

f or (ctr2=0 ; ctr2<no_neurons [ctrl] ; ctr2++) i 
free(taps[ctrl] [ctr2] ) ; 

} 

} 

for  (ctrl=0 ;  ctrl<no_layers  ;ctr  !++)•( 
f ree(taps[ctrl] ) ; 

> 

free(taps) ; 

return ; 

} 
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